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SUMMARY 


By employing well-known relations connecting conjugate har- 
monic functions, integral equations are obtained for the incom- 
pressible velocity distribution over thin airfoils, isolated and in 
cascade. Solution of these integral equations with the boundary 
condition v = u (dy/dx) in place of the usual v = V (dy/dx) |[i.e., 
classical] condition gives distributions free of the usual leading- 
edge singularity. When the classical boundary condition is 
introduced, the integral equations reduce to the familiar isolated 


thin airfoil theory and its extensions to cascade. 


(1) INTRODUCTION 


bw CLASSICAL THEORY of potential flow past iso- 
lated thin airfoils has proved a valuable tool for 
the analysis and prediction of the aerodynamic char- 
acteristics of wings and airfoils. Originally introduced 
by Birnbaum, Munk, and Glauert in several different 
forms, it has been developed and refined by many 
later workers. Recently, renewed interest has been 
stimulated by the study of higher order theories. 
Lighthill' and Van Dyke*® have shown that very accu- 
rate results can be obtained by simple modifications of 
the classical first-order thin airfoil theory. 

Thin airfoil theory for cascades is as fundamental in 
the aerodynamics of axial-flow turbomachinery as is 
Birnbaum-Munk-Glauert theory for the aerodynam- 
There is a very extensive literature on 
In particular, 

One is the 


ics of wings. 
potential flow through airfoil cascades. 

two approaches have been widely used. 
singularity-distribution method, in which sources and 
vortices are distributed to represent each airfoil. The 
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integral equations which result are solved by expressing 
the singularity distributions in the form of series and 
satisfying the boundary conditions at a finite number 
of points, thus reducing the problem to the solution of 
a number of simultaneous equations for the coefficients 
of the singularity distribution. This approach has 
been used by Pistolesi® and Lieblein‘ among others, 
and has recently been discussed extensively by Scholz* 
and Schlichting.® 

The second approach is that of conformal transfor 
mation, in which the methods of complex variable 
and conformal transformation theory are applied to 
the boundary value problem represented by the poten- 
tial flow through a cascade of airfoils. This approach 
has been used for calculating the exact flow about air- 
foils of special shapes as well as the exact and approxi- 
mate flows about arbitrary airfoils. A review of several 
such methods is given in a report by Tyler.’ In par- 
ticular, Garrick* has shown how to find the exact solu- 
tion for the flow past a cascade of arbitrary airfoils by 
extending the Theodorsen-Garrick method for isolated 
airfoils to cascades. However, his solution is an iter- 
ative one which requires extensive numerical work. 

The present paper describes a method of finding flow 
past arbitrary airfoils of small (but finite) thickness and 
camber either isolated or in cascade which tries to strike 
a useful compromise between accuracy and numerical 
labor. It is based on the complex variable approach. 
As such, it has many points of similarity with the classi- 
cal Munk theory for isolated thin airfoils. Like the 
Munk theory, it leads to explicit formulas for lift 
moment, and velocity distribution in terms of the air- 
foil shape when the boundary condition of classical 
thin airfoil theory is used. In the case of unstaggered 
cascades, such formulas have recently been published 
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Fic. 2.1. Isolated airfoil with chord and circle planes. 


by Dorr® and Nickel,'’ but for staggered cascades the 
results derived here appear to be new. 

These Munk-type integral formulas also form the 
basis, hitherto unavailable, for the kind of higher order 
refinements which Lighthill' and Van Dyke’ have de- 
veloped for isolated airfoils. It is hoped that these 
possibilities will be explored in the near future. 
Application of a boundary condition somewhat dif- 
ferent from the classical one leads to an integral equa- 
tion for the velocity distribution, which, when solved 
by iteration, yields results which compare very well 
with exact solutions. In particular, the very large 
velocities in the neighborhood of the leading edge, pre- 
dicted when the classical thin airfoil boundary condi- 
tion is used, are avoided. Furthermore, the agreement 
with exact solutions is as satisfactory near the leading 
edge as it iselsewhere. Below the stall, agreement with 
experiment is also quite good. In view of its success 
in predicting aerodynamic characteristics of cascades, 
the present theory is of value in the continuing effort 
to improve the efficiency of turbomachinery. For a 
given airfoil, the solution of one singular integral equa- 
tion will give the velocity distribution at all angles of 
attack with good accuracy, even near the leading edge. 


(2) MrtTHoD OF ANALYSIS 


The first objective in the method to be presented 
here is to obtain a pair of integral equations which re- 
late on the airfoil contour the Cartesian velocity com- 
ponents u and v of the two-dimensional incompressible 
potential flow past the airfoil section. These relations 
are consequences of the fact that u — iv is an analytic 
function of the complex variable z in the region exterior 
to the section. They are basic to the solution of both 
the direct (given an airfoil, find the velocity distribu- 
tion) and the indirect (given a velocity distribution. 
find an airfoil which produces it) problems of thin 
airfoil theory. The method applies equally well to the 
isolated airfoil and to the staggered cascade of airfoils; 
for convenience, however, the discussion will be re- 
stricted at first to the former. 

The integral equations which are being sought are 
in their simplest form when expressed for the region 
exterior to a circle. It is for this reason that a trans- 
formation is made from the airfoil plane to a circle 
plane (the ¢ plane). It is desirable, however, to avoid 
the complexities arising in the construction of the con- 
formal map of the airfoil exterior onto the exterior of 


the circle. For this reason the usual approximation of 
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thin airfoil theory is introduced: _ the airfoil is replaced 
by a straight line segment (its chord) and the velocity 
components u and v at the surface of the airfoil ate 
transferred to this line segment (see Fig. 2.la). It js 
the complex plane, slit along this chord (see Fig. 2.1b) 
which is now mapped onto the exterior of the circle 
¢| = a (Fig. 2.1c) by the Joukowski transformation, 

The velocity components of the transformed potential 
flow in the ¢-plane are denoted by u* and v*.  Appli. 
cation of the integral relationships directly to «* and 
v*, however, does not give the desired form when ex. 


Instead, we consider the rea] 
* 


pressed in terms uw and v. 
and imaginary parts of the product — i*)¢, 
function whose physical significance lies in the fact that 
its integral around the circle is proportional to the 
circulation around, and net outflow of fluid from, the 
circle. It should be noted in passing that the chord 
in the z plane (see Fig. 2.1b) and the circle in the ¢ plane 
are not streamlines of the flow. 

The function (w* — iv*)¢ is continuous on the circle 
¢| = aand analytic outside it except for a pole at ¢ = 
o, If the principal part, (w..* — iW..*)f, at this pole 
is subtracted from the function, the difference is ana- 
lytic and bounded in ¢ >a. The integral equations 
then result from an application, to the function /(¢) = 
[(u* — u.*) — i(v* — v.*)]¢, of the following theo- 
FUNDAMENTAL PRINCIPLE. /f f(¢) ts analytic 
¢| = a, continuous on this circle 


rems: 
exterior to the circle 
and bounded as § > @©, then 


= (1/2mi) { flae’’) cot — @) + C 
0 
(2.1) 
where C= (1 an) flae’®) do = f(@) (2.2) 
0 


This principle, which is associated with the names of 
Hilbert, Dini, Fatou, and others, is derived in the Ap- 
pendix by elementary methods. The companion re 
sult for the interior of the circle, also proved in the 
Appendix, is given as: /f f(¢) is an- 
alytic in the interior of the circle |¢| = a and continuous 
on this circle, then 


flae®) = —(1/2zmi) flae’®) cot [(@ — 
0 
(2.3) 
where C= (1 an) f(ae’*)dp = f(0) (2.4) 
0 


The two integral relations thus obtained are actually 
two different forms of the statement that « and —? 
are conjugate functions. They lead, respectively, to 
the solutions of the direct and indirect problems when 
they are supplemented by the boundary condition of 
potential flow at the airfoil surface, which relates 
u and v in terms of the slope of the airfoil profile. 

It should be pointed out that the present paper 1s 
concerned primarily with the direct problem of a given 
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Consideration first of 


airfoil, isolated and in cascade. 
the case of the isolated airfoil demonstrates the appli- 
cation and power of the method without obscuring them 
Here the familiar results 


in algebraic computation. 
for the direct problem appear as special cases and new 
insight is provided into the implications of the classical 
solution. The method is then applied to a staggered 
cascade of thin airfoils, leading to new results. 


(3) THe ISOLATED AIRFOIL 


An isolated thin airfoil with small camber and chord, 
rc. as shown in Fig. 2.1, is placed in a stream of velocity, 
l’, at angle of attack, a. The problem of determining 
the flow about the airfoil is a boundary value problem 
of potential theory. In accordance with the thin air- 
foil approximation, the boundary-value problem is trans- 
ferred to the chord line —c 2 Sx Sc 2,y = 0. On 
that line, the Cartesian velocity components u and v 
satisfy the continuity and irrotationality conditions 
and also the boundary condition that the flow be tan- 
gential to the airfoil surfaces. 

It is well known that for two-dimensional incom- 
pressible flow, methods of complex variable theory may 
The continuity and irrotationality conditions 
iv is an 


be used. 
imply that the complex velocity w = u — 
analytic function of the complex position variable 
s=x+/y. The problem is to find a relation between 
wand v, on the chord line —c 2 Ss: 0, which 
results from this fact. A second relation between u and 
’, expressing the tangential flow condition, is avail- 
able. Combination of these two relations will permit 
us to eliminate either the w or the v component, ex- 
pressing the other one entirely in terms of the airfoil 


shape. 


Transformation to the Circle Plane 


Assume initially that the airfoil contour (see Fig. 
2.1a) has been mapped onto the chord line —c 2 S$ 
x Sc 2, y = 0 (see Fig. 2.1b) by the analytic function 
which transforms the plane exterior to the contour 
conformally into the plane slit along the chord, trailing- 
edge points being made to correspond. It is this 
transformation from the slit plane to the contour plane 
which is subsequently approximated by a perpendicular 
projection. The value of « — iv at a point in the slit 
plane is taken equal to its value at the corresponding 
point in the contour plane. This is in contrast to the 
usual procedure since the velocity u — iv is not being 
modified under the conformal map. As a result, the 
slit is not part of a streamline, but the function « — iv 
is an analytical function of z. This property serves to 
establish a relationship between u and 2. 

Next map the region exterior to the chord conform- 
ally onto the region exterior to the circle ¢ = a (see 
Fig. 2.1b) by the Joukowski transformation 

(3.1) 
On the circle itself, ¢ = ae”, so that 


2a cos 6, 


x= 


THIN 
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The parameter a is thus to be identified with the 
quarter-chord, a = c,4. The derivative of the trans- 


formation is 
= 1 


which has the value 


- 
dz/dt = 2ie~” sin = ie~" 
on the circle. The notation g(@) is introduced in place 
of sin @ to indicate an analogy with the cascade case 


which will appear later. 


The Cartesian velocity components u*, v* in the ¢ 
plane are given by 
u* — w* = (u — iv) (dZ de) (3.5) 


where, as indicated above, « — v/v is the velocity in the 
slit plane. On the circle, using Eq. (3.4), we can then 
write 


(u* — iv*)t = (u* — iv*)ae” = 
2ai (u — iv) (3.6) 


Applying the Fundamental Principle on the circle 


¢ = ae” to the function ¢[(u* — u..*) — i(v* — iv.)); 
that is, to ¢(u* — w* — Ve‘); gives 


ae” [u*(0) — iv*(0) — Ve '*] = (1 2m) 


[u*(o) — iv*(o) — Ve X 
cot [((@ — ¢) 2)|}d@+C 


The Fundamental Principle Corollary, Eq. (2.5), 


plied to f(¢)=¢ = ae", yields 


4 


(1 ae’® cot [(@ — @) 2]d@ = —ae" 
0 
since /(¢) = Oat ¢ = 0. Using this and Eq. (3.6) in 
Eq. (3.7), we find 


2ag(@) [v + iu] — = (1 2ni) | 2ag(d) X 


+ tu(d)] cot ((@ — 


0 


(33.8) 


The constant C can be related to the total circulation 
I’ (positive clockwise) and net source strength Q of the 


flow. These are given by 


—r+iQ=4+ (udx + +i f (—vdx + udy) = 
S (u — = f (u* — iv*)de = 


aie’ [u*(0) — (3.9 
0 


(which also expresses the fact that T and Q are the same 
in Z plane and ¢ plane). Now Eq. (3.7) is integrated 
with respect to @ from 0 to 27. Since the Cauchy 
principal-value of the integral of cot(@ — @) 2 has the 
value zero, we find upon changing the order of integra- 
tion that 


2rC = ae’ [u*(@) — iv*(0) \dé 
0 


and, consequently, 


— 
) 
| 
X 
(3.7) 
») 
(2.1) 
(2 | 
«C4 
| 
»+C 
( » 3) 
(2.4) 
ually 
1 
vy, to 
am = 0 2) 


JOURNAL OF THE AERO 


Yt 

Yt 


X CHORD 


Fic. 3.1. Section of airfoil with coordinate system. 


(3.10) 


If Eq. (3.10) is substituted into Eq. (3.8) and the real 
and imaginary parts are separated, the result is 


g(@)v(0) = I’ cos (@ — a) +(1 2am) X 
( 


) 


u(d) cot — ¢) + (Q (3.11a) 


= I sin (@— a) — (1 2am) X 
0 


v(m) cot — @) 2) d@ + 42a) (3.1 1b) 


These two equations express the relation between 1(@) 
and v(@) on the chord line, where position, x. is related 
to 6 by x = 2a cos @ witha = c 4. The constant Q is 
determined by the condition that there be no flow across 
the contour. Thus Q = 0. The circulation is found 
from the Kutta condition that the velocity be finite at 
the trailing edge, @ = 0. Since g(@) = O, this condition 
can be satisfied only if the right-hand side of Eq. (3.11b) 
vanishes at 6 = 0. The circulation is thus found to be 
replacing g(@) by sin 0, a by c/4| 


V sin a — (1/27) X 
sin v(@) cot (¢ 2) de | (3.12a) 
0 


Alternatively, integration over (0, 27) in Eq. (3.1 1b) 
gives 


= u(@) sin 6 (3.12b) 

0 
The velocity relations may therefore be written 
2(@) = cos (@ — a) sin 6] + (1/27 sin X 

sin u(@) cot [((@ — ¢) 2] dd (3.13a) 

0 
u(@) = — cos sin 6] sin a + cos a} 
(1, 27 sin 6) sin @ v(d) X 
0 


[(@ — $)/2] + cot (¢/2)} dé (3.13b) 
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Eqs. (3.13) are the final expressions for the relations , 


between u and expressed in the circle-plane coord. 
nate 
plane coordinate x can be found by using v = (¢ 2 


The corresponding expressions in the physical. 


cos 6. 

Eq. (3.13a) is useful for solution of the so-called jn- 
direct problem in which the airfoil shape is to be found 
from a given distribution of velocity wu. 
12 for a discussion of the indirect problem. ) 


(See reference 
The pres- 
ent paper is concerned with the direct problem in which 
the airfoil is given and the velocity is to be found, 
Here, Eqs. (3.12) and (3.13b) are of primary interest, 
In order to complete the solution of the direct problem 
v must be eliminated from these equations by use oj 
the boundary condition that the flow be tangent to 
the surface of the airfoil. 


Thin Airfoil Approximation and the Classical Solution 


It is important to notice that the Eqs. (3.1 la) and 


(3.11b) can be considered to be exact. The “thin air- 
foil’ approximation is made as indicated in Section (3 
by assuming that the mapping of the airfoil section 
onto its chord line is simply a perpendicular projection 
onto this line. If y = } (x) is the double-valued func. 
tion describing the airfoil contour, the exact boundary 
condition is 


v= u(d¥/dx) onthe airfoil contour = (3.14 


This is to be adapted to the thin-airfoil approximation 
by the indicated transfer to the chord. 

In the classical solution of the direct problem, the 
further simplification is made that “ may be replaced 
in Eq. (3.14) by I’ cos a, the horizontal component 
of the undisturbed flow. Thus, the classical boundary 
condition is 


v = Vcosa(dY/dx) ony = 0 (3.15 


Actually, cos @ is usually replaced by unity since a is 
small, but it is convenient to retain the term here to 
show explicitly the sinusoidal dependence of the solu- 
tion upon a. 

We shall now show how the equations of classical 
thin airfoil theory follow from the present approach 
and the use of the boundary condition (3.15). First 
it is convenient to express the velocity components u and 
v as perturbations of the free-stream flow: 


u= Veosatu’, v= Vsinata’ (3.16 


Egs. (3.12) and (3.13) become 


v'(0) = (1/27 sin 6) u’(@) sin 
¢ cot — 2] (3.1%) 
u'(@) = (—1 27 sin @) v’(@) sin 
cot — ¢) 2)d¢ + (1 cx sing) (3.17b 
r= —(c 2) v'(d) (1 + cos dp (3.15a 
= (c/2) u'(p) sin ¢ (3.18b) 


re 
re 
ul 
Yu Cc Ye 
} 
v 
Si 
P 
is 
u 
| 
a 
1 


relations , 


€ coord)- 
physical. 


= (¢? 


‘alled in- 
be found 
“eference 
‘he pres- 
in which 
> found, 
interest. 
»roblem 
Vy use of 
igent to 


ition 


la) and 
hin air- 
tion (3 
section 
Djection 
fune- 
yundary 


(3.14 
mation 
‘in, the 
eplaced 


ponent 
undary 


(3.15 
ce ais 
1ere to 
e solu- 
assical 
proach 

First | 
smand } 


(3.10 


3. 17a) 


[Na 


8b) 


A THEORY OF THIN AIRFOILS 413 
Let the upper and lower surfaces of the airfoil be thickness of the airfoil. The upper surface of the air- 
represented by the functions y = y,(x) and y = y,(x) foil, which corresponds to the upper side of the chord, 


respectively. Thus, (see Fig. 3.1) is represented by the upper half of the circle, 0 S 6S rz. 


The lower surface is similarly described by 7 S @ S 2r. 
Yu = MM The integrals in Eqs. (3.17) and (3.18) may now be 
separated into integrals over these intervals with the 
where y = Y-() represents the camber line and 2y, the corresponding boundary condition [c.f. Eq. (3.15)] 
introduced. 
u'(6) = —(1 27 sin 6) sin a + cos a(dy, dx)| sin ¢ cot [((@ — ¢) 2] + 


9 


[I sin a + I’ cos (dy,/dx)] sin ¢ cot — ¢) 2] doy = [(1 — cos @) sin 4] sina — cosa mr) X 


0 


(dy, dx) [1 + cos @) (cos @ — cos doy cosa m) | (dy, dx) [sin @ (cos — cos 6)| de 
0 


Thus u’ may be separated into odd and even components, u,’ and u,’, where u,’ represents the effect of thickness and 
u,/ represents the effects of camber and angle of attack. 


u,/(@) = [(1 — cos @) sin 6] sina — (I’ cosa (dy, dx) [((1 + cos @) (cos — cos 9) doe (3.19a) 


0 
u,/(0) = —(V' cos a@ m) { (dy, dx) [sin @ (cos @ — cos 6)| do (3.19b) 
Similarly, 
= (1 27 sin @) + u,'(d)] sin cot — ¢) 2] de — 
[—u,.'(@) + u,'(o)] sin cot — ¢), 2] dot = (1/7) X 


[sin @ (cos — cos + (1 sin 8) u,'(@) [sin® @ (cos — cos 4) 


0 
so that = (1/m) lsin @ (cos — cos 6)| do (3.20a) 
0 
is even and v,/(0) = (1 6) u,'(@) @ (cos — cos de (3.20b) 
0 
isodd. The circulation becomes (T, = 0) 


= sin a — (I’ cos a m) (dy, dx) (1 + cos ¢) ia | = u,’ sin @ dd 
0 0 


Changing to rectangular coordinates in the chord plane x = (¢ 2) cos @, & = (c 2) cos ¢, gives for camber 


v-'(x) = (1 u,.(é) [dé —x)], = 2f u,'(&) dé 


u(x) = 2) — x} /[(e/2) + x] I’sin a — cosa X (3.21) 
(dy, dx) Vi(c 2) + [(e 2) — [dé 
and for thickness 
= V(c2/4) — u,'(&) V (ce? 4) — & [dé 
(3.22) 


= —(V cosa 7m) (dy, dx) [dé — x)] 

These are the familiar relations of classical thin airfoil 
theory. The camber results are given, for example, in A single set of equations leads directly to the relations 
reference 11, pp. 45-47, for a = 0 and both sets are for thickness and camber effects in the direct problem, 
giving solutions to the integral equations as fomulated 


given in reference 13. 
This point of view further provides 


The compactness of the present method is apparent. by Soéhngen.'* 


Fic. 4.1. Cascade of thin airfoils. 


insight into the significance of the separation of camber 
and thickness effects which will be indicated in the next 
section. 


The Boundary Condition 


It may be noted from Eq. (3.19) that in general, 
u’(0) becomes infinite at the leading edge, @ = 7. This 
is characteristic of classical thin airfoil theory and 
makes it unsatisfactory for use in finding pressure dis- 
tributions near that point. In addition, it should be 
remarked that in the classical solution for an airfoil, in 
which both thickness and camber are considered by the 
simultaneous distributions of sources and vortices along 
the chord, the source strength is determined from thick- 
ness alone and the interaction effect upon it of camber 
is ignored. 

These considerations led C. B. Smith in 1949 to the 
introduction of the concept of ‘‘one-sided”’ source (c.f. 
reference 15). This device permitted the interaction 
of thickness and camber and eliminated the leading- 
edge singularity when used in conjunction with the 
boundary condition (3.14), applied at the chord, in- 
stead of Eq. (3.15). The present method is largely a 
result of investigations based on this concept. 

If the boundary condition 


v = u(dy dx) on the chord (3.23) 


is introduced into Eq. (3.13b) a singular integral equa- 
tion in results: 
u(0) = V}[(1 — cos @) sin 6] sin a + cos al X 
(—1/27 sin @) u(@) (dy dx) sin 
¢jcot [(@ — 2] + cot (¢ 2)'dp (3.24) 


where dy dx is evaluated at the x associated with ¢ by 
x = (c/2) cos @ Extensive solutions of this integral 
equation have been carried out employing finite-differ- 
ence methods on a digital computer both for the isolated 
airfoil under discussion and for cascades. The numeri- 
cal solution (discussed later) is found to converge very 
rapidly, three or four iterations usually being sufficient. 
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Comparison of boundary conditions is of interest 
here. Writing Eq. (3.14) in terms of perturbation ye 
locities and camber thickness components gives 


V sina + 9,’ + = (V cosa + 
ue’ + uz’) dx) + (dy, dx)] 


where the odd function y, is equal to y, on the upper 
surface, —y, on the lower. These give, upon sepa- 
ration into odd and even parts, 


(V cosa + u,’) X 

(dy. dx) + u,.'(dy, dx) 
(lV cos a + u,') X 
(dv, dx) + u,'(dy, dx) 


sin a + 2%, 


« 


in contrast to the relations used above in the classical 
solution, 


I’ sin a + 2,’ = cos a(dy. dx). 


v,’ = V cos a(dy,/dx) (3.26) 
Finally, it may be shown, by the use of Fourier Series, 
for example that the integral equation (3.24) has a 
continuous solution “(@) and that this continuous solu- 
tion is unique. There is no singularity in “(@) at the 
leading edge and the pressure peak immediately behind 
this point is seen in the plots presented subsequently. 
There is another important consequence. The function 
— (0° Oa’) is also a solution of Eq. (3.24), and, by the 
uniqueness property, 


(O°u Oa?) + u = O 
Thus, = cos a + sin (3.27 


and it is necessary to solve Eq. (3.24) for only two 
values of @ in order to get the solution for arbitrary a 


(4) CASCADE OF THIN AIRFOILS 


The use of the method of conjugate functions in thin 
airfoil theory has been examined above in some detail 
for the isolated airfoil. The extension to cascades of 
thin airfoils, although somewhat more involved alge- 
braically, is immediate and follows the same steps. 
The resulting relationships can be reduced, as a special 
case, to the extensions to cascades made by H. Schlicht- 
ing® for the classical theory. 

A cascade of thin airfoils with small camber is shown 
in Fig. 4.1. The airfoils have chord c, gap 7, and are 
staggered at angle 8. The boundary value problem is 
again transferred to the chord of each blade, and a 
conformal transformation is used to map the plane, 
slit along these chords, onto the exterior of a circle 
An application of the Fundamental Principle produces 
the relations between « and v. The only complication 
in the cascade case (other than the algebraic one 
arises because of the presence of singular points in the 
circle-plane outside the circle. However, it is found 
that these can be handled by a consideration of the 
velocities infinitely far ahead of and behind the cascade. 
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Transformation to the Circle Plane 


A staggered cascade of straight line segments (see 
Fig. 4.2a) can be transformed into a circle (see Fig. 
| 2b) by means of the mapping function given in refer- 
ence 11, p. 95. 


The logarithmic nature of the transformation makes 2 


a multiple-valued function of ¢ so that the circle. 
- = ae’, is transformed into the infinite set of line seg- 


ments 


n=0, £1, +2... (4.2) 
where Y denotes position of the line segment which lies 
along the x axis. Also note that, for 8 = 0, if k > Oas 
7— © in such a way that kr wa > 1, then Eq. (4.1) 
reduces in the limit to the isolated airfoil transforma- 
tion (3.1). The real parameter k (0 < k < 1) has no 
simple physical interpretation and will be determined 
later in terms of 8and 7 c. 

On the circle ¢ = ae’, the relation between 2 and @ 
is given by Eq. (4.2), where the function X (@) is 


X(6) = (7 7) [cos @ arc tanh (cos 6A) + 


sin 8 arc tan (sin 6 A*)] (4.3) 


Here the real parameters A and A™* are related by k by 


K=(1/k + &)/2. K* = (1/k — = 


(A (4.4) 
The derivative of the transformation 
dt (1 k) + (¢ a) 
|+ 
(1k) — (fa) (€/a) +k 
] 
| (4.5) 
(C/a) —k f 


which becomes, on the circle, 


Z- PLANE C- PLANE 
(a) (b) 


FiG. 4.2. Transformation of staggered cascade of straight line 
segments onto a circle. 


THIN 
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dz dt = dx tae” dé = 2ie~"G(0) 


where 
(4.6) 
r A siné@cos sg — A* cos @sin 
G(@) = 


2ra K? — @ 


Now the leading and trailing edge of the straight line 
segments correspond to the points ae’ and ae“” on 
the circle where dz dé = 0; 6; is the angle of smallest 
absolute value satisfying [see Eq. (4.6) } 


tan 67 = (A* A) taneB (4.7a) 


and 6; = 6, + m The upper side of the line —c 2 S 
x < c 2 corresponds to 67 < 6 = 6,, while the lower 
part corresponds to 6; S 6 S @7 + 27. 

A useful set of relations between 67 and 8 can be ob- 


tained from Eq. (4.7a) 
sin = K* sin B/J, J = 
(4.7b) 


cos 67 = K cos 8 J, 

V K? — sin? 8, — cos? = (KK* 
With these, Eq. (4.6) can be put into the more compact 
form 


— (dx /2ad0) = G(6) = (Jr X 


[sin (@ — 67) (K? — cos? 6)] (4.8) 


(Or) X(0,). 
Then from Eqs. 


The length of the chord is ¢ = 
Note that X(6 + 7) = —7X(@). 
(4.3) and (4.7) it follows that 


c tr = (2 m) [cos 8B arc tanh (cos 67 A) + 
sin 6 are tan (sin 67 A*) | 
(2 m) [cos 6 are tanh (cos 8 J) + | 
sin 8 are tan (sin 8 J/)| 


(4.9) 


This equation relates the transformation parameter 
K (and thus k) to the two physical parameters ¢ 7 and 
3. Tabulations of this transcendental relation are 
available in the literature, in somewhat different nota- 
tion. (See, for example, reference 8, pp. 3 and 4, where 
k-! = b/a, K = cosh y, K* = sinh y.) 


Application of the Fundamental Principle 


The Cartesian velocity components in the circle 
plane are related to those in the cascade plane by 


u* — iv* = (u — iv) (dz/dg) (4.10) 
On the circie, then, from Eq. (4.6), 
t(u* — iv*) = ae” (u* — iv*) = 
Zia (uw — iw) (4.11) 
Here ¢(u* — iv*) is bounded as{— ©. From Eq. 
(4.1) it is seen that the point ¢ = © corresponds to the 


points s = (1 2)ir(2n + 1)e~ where n is any positive 
or negative integer, or zero. These points all lie on the 
center line of the cascade, halfway between adjacent 
chord lines, and the velocity « — @v is finite there. 
Further, note from Eq. (4.5) that ¢(dz df) — 0 as 


¢—~o. Therefore, 


lim ¢(u* — iv*) = lim (u — iv)¢(dz dg) = 0 (4.12) 


t—« 
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(3.27) | 
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Fic. 4.3. Velocity triangles for the cascade. 


The function ¢(“* — 1i*) is not, however, analytic 
everywhere outside the circle lz | = a. Reference to 
Eq. (4.5) indicates that dz df has simple poles at the 
points ¢ = +a k. Substitution of these values into 
Eq. (4.1) shows that they correspond to points infi- 
nitely far in front (v = — ©) and in back (x = +) 
of the cascade, respectively. At the former point let 
the complex velocity be Ve~', at the latter l.e7'*.2 
(It is well known that a cascade turns the flow so that 
the velocities far in front and far in back are not the 
same. Notice that the angles are defined with respect 
to the chord.) Then, from Eqs. (4.10) and (4.5), 
follow the relations 


+ (€/a) 
(u* — iw") => 


2nra — (¢/a)’ 

In view of these relations it is seen that the function 

27a 


k) + (¢/a) (1/k) — (¢ 
is both analytic outside the circle and bounded as ¢ > 
and the Fundamental Principle may be applied to the 
function F(¢). 
For convenience, define a function //(¢) to be 


Vee~*™ | 
(1/k) — (¢/a) 


(4.14) 


Me) re | Ve ™ 
+ (€/a) 
so that F(¢) is given by 
= g[u* — w* — 
From the Fundamental Principle, Eq. (2.1), it follows 
that 


ae’’ [u*(@) — iv*(0)] — ae“ Hae") = 
(1/271) ae’®[u*(¢) — iv*(o) — 
0 
H(ae’®)| cot — ¢)/2]deé + C 


(4.15) 
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The function ¢ //(¢) is regular inside the circle anq\ 
vanishes at the origin. 

Therefore, the Fundamental Principle Corollary, Eq. 
(2.5) can be applied to yield 


(1 (272) cot [(@ — 2]d¢ = 
0 
—ae™ (ae" 


Using this and Eq. (4.11) in Eq. (4.15) gives, in analogy 
to Eq. (3.7), 


2aG(0) [v(@) + iu(@)] — 2ae” Hae”) = 
(1 ani) 2aG(o) [v(@) + x 
0 


cot [(@ — 2]d@ + C (4.16) 


Simplification of the Velocity Relations 


The definition of total circulation T° and net source 
strength Q given in Eq. (3.9) holds for an airfoil in cas- 
cade also. Therefore, if Eq. (4.15) is integrated on 
6 from 0 to 27, the result analogous to Eq. (3.10) is 
obtained for the constant C: 


= 7 2ae” H(ae")da 

0 
The integral on the right is easily evaluated by letting 
¢ = ae", dé = aic'd9, whence it becomes a contour in- 
tegral, around the circle, of —2///(¢). Since //(£) is 
singular only outside the circle, the integral vanishes 
by Cauchy's theorem, and it follows that 


C= (0+ (4.17 


asin Eq. (3.10). 

Another expression for C can be obtained from Eq. 
(2.2), where it is shown that C = F(). Letting 
¢— o in Eq. (4.13) and noting Eq. (4.12) yields 


C = —(re~ (Vie ** Voe~*™*) (4.18 
When this expression for C is separated into real and 
imaginary parts and equated to the right-hand side 


of Eq. (4.17), it is seen that 


Q = 


cos (ag + 8) — cos (ay + 8)]  (4.19a) 


sin (ag + B)] (4.19b 


sill (ay + B) — Vo 
Reference to Fig. 4.3 shows that Q 7 is the difference 


between the components of I; and I) normal to the 
If flow continuity is to be satisfied, 
The quan- 


cascade line. 
this must of course be zero so that VY = 0. 
tity I r is the difference between the components of 
I’, and I’, parallel to the cascade line and represents the 
turning of the flow through the cascade. An expres- 
sion for the turning angle in terms of T is found from 
Fig. 4.3. First it is convenient to introduce the vector 
mean I’e~‘* of the velocities far upstream and down- 
stream of the cascade: 


= (1/2) + (4.20) 


This velocity vector is shown in Fig. 4.5. In terms ol 
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it. the turning angle © is found, with the help of Fig. 
43 to be 
— ao = arc tan Cr X 


'2 cos (a + 8) [1 — (I 2v7)?]}) 


G0 = 


(4.21) 


From 0 = 0 it follows that C = iI 27 so that by 
Eqs. (4.18) and (4.20) 


Ve C= Ve — 
Thus 
+B) » 210 
Hae”) = (vr mak) [(1/k2) — + 


Eq. (4.16) can now be written in final form. Using 
Hae) from Eq. (4.22) and the fact that C = iT 27, 


separating into real and imaginary parts gives 


= G(o)u(@) cot do — 
2e Jo 2 
sin 26 


Sma K? — cos? 6 27a 


AK* cos 6 cos (a + 8) + A sin @ sin (a + 8 


(4.25a) 


A? — cos? @ 
wi 
G(6) u(@) = — cot dg + 
2r Jo 2 
KK* Vr 


tra K* — 0 
K sin 6 cos (a + 8) — A* cos 6 sin (a + B) 
K? — cos? 6 


(4.23b) 


These equations are the velocity relations we are 
seeking, analogous to Eqs. (3.11) of the isolated case. 
Eq. (4.23a) involves an integral over u, and as such 
is useful in the indirect problem where the velocity 
distribution is given. In the direct problem being 
considered here, where the airfoil is given, it is Eq. 
(4.23b) which is of primary interest. The circulation 
lis determined from the Kutta condition of finite ve- 
locity at the trailing edge, 6 = 67. Since Eq. (4.8) 
shows that G(@r) = O, the right-hand side of Eq. 


(4.23b) must also vanish for @ = 47. Thus, 
u() | jsin (@ — | sin + sin @ kK? 
| Jo (@ — Leos — cos 
dy 
V dx 


The remarks at the end of Section (3) concerning the 
integral equation (3.24) apply equally well here: Eq. 
(4.28) has a continuous solution so that the pressure 
peak immediately behind the leading edge may be 
This solution is again 


used in design procedures. 


THIN 


> — cos? @) — 


sin a + AJ sin 6 cos (a + B) — A*J cos 6 sin (a + 
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r= (2Vrsin a/J) + (2KK*/ J*) X 


f aG(@)v(@) cot [(@r — @) 2) d@ (4.24) 
0 

Write «# and v as perturbations of the mean flow 
Ve""*: u = Veosat+u’,v = Vsna+v’. Then 
Eqs. (4.23) become simply 


G(0)v'(0) = (1 2m) cot — de — 
0 


(T Sara) [sin 26 (AK? — cos? 6)]  (4.25a) 
G(@) = —(1 2x) | G(@)u'(@) cot 
0 
— 2] + (1 42a) X 
[AA* (AK? — cos? @)]  (4.25b) 
where 
(ER rf aG(@)v'(@) cot 
0 
— o) 2) d@ (4.26) 


An expression for the moment about the airfoil center, 
positive clockwise, can be derived from «’(@) by integra 
tion. The pressure disturbance in linearized theory is 
related to u(@) by p = 


where p is the fluid density. 


VW = (Prop Povottom) xdx 


f top — bottom) Xdx 


0 


—pl(u — Vcosa) = —plu’ 


Thus, the moment is 


ll 


since dy = —2aG(@)dé from Eq. (4.8). Substituting 


for u’(@) from Eq. (4.25b) gives 


M pV = (1m) X 
0 


[sin 8 are tanh (cos ¢ A) — cos 8 are tan 


(sin @ A*)|]d@ (4.27) 


Eqs. (4.25b), (4.26), and (4.27) are the velocity rela 
tions for the direct problem of thin airfoils in staggered 
cascades. It is necessary only to relate v to the airfoil 
profile by means of the boundary condition of tangential 
flow in order to complete the solution. With the 
modified boundary condition Eq. (3.23), Eqs. (4.25b) 
and (4.24) give the integral equation 


sin 67 + sin @ 


(K? — cos? ax) |t 
cos @ — cos Or f 


(4.28) 
J* sin (0 — Or) 


unique and has a sinusoidal dependence upon the angle 
For v = 0 (i.e., for v’ = V sin a), these equations 


become the corresponding relations for the flat plate 


x 
| 
—_ 
| 
| 
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cascade (reference 10, pp. 95 ff.) : and (3.13). The moment also reduces to the know, — T 
‘ — isolated airfoil expression. 
trp = Vceosa+ (V sin a/J) X I 
AA* tan [(9 — Or) 2] — cos B sin The Classical Boundary Condition 
Typ = sina J If, instead of Eq. (3.23), the boundary conditio; 
Mep pV = (2V wr? sin a ln (K K*) (3.15) of classical thin airfoil theory is inserted int 
Eqs. (4.25b), (4.26), and (4.27), explicit expressions 
where the subscript FP refers to ‘‘flat plate.” foru. Tl, \/ are obtained. Thus. 
As7— with 3 = 0, the relations (4.25) and (4.26) : 
reduce to the analogous ones in the isolated case, (5.12) v= Vsina +0 a 
I’cosa (dy dx) on y=0 (4.29 
gives 
2s 
— = cosa an) f (dy dx) [G(@) G(9)] cot [(@ — 2] de + 
0 
(KA* 27J) pp sin (6 — 67) | 
2r 4.5 Ja 
= cosa an) (dy dx) [G(¢) G(@)] {cot [(@ — ¢) 2]de E 
0 
[(K? — cos? 67) (K?2 — cos? 6)] — do 
Ce 
—Trp = — (tKA*V sin a (dy dx) | [1 + cos (@ — 6r)] (K? — cos? de (4.30b 
0 
» SO 
M — Mep pi = (2avr cos a m) G(@) (dy dx) [sin 8 are tanh (cos A) — cos Bare tan (sing A*)|dd  (4.30¢ 
0 
Here x and @ (or @¢) are related by x = (8) where - —- 
X (0) is given by Eq. (4.3). Eqs. (4.30) represent an both the top (@r < < + =) and the bottom (4; + 
explicit solution of the direct problem for thin airfoils r <6 < 67 + 2m) so that airfoils of arbitrary shap: 
of small camber in staggered cascades. As such, they can be handled (subject to the small thickness and cam- | 
are an exact analogy to the Munk integral formulas of ber restriction, of course). 
isolated thin airfoil theory. In particular, Eq. (4.50a) 
is the “inversion formula’ corresponding to the Sdéhn- Reduction to Cartesian Coordinates pal 
gen formula for an isolated — Comparison be- For nonzero stagger, separation of Eqs. (4.25) int 
siaatas results based upon Eq. (4.28) and those obtained odd and even components does not simplify the funda- | 
from the classical sone of solution [Eq. (4.30) ] will be mental relations as it did in the case of the isolated air | 
given in the last section. oom foil. Here the transformation , = X(6) of Eq. (43 Fe 
It should be noted that here “(@) will, in general, be is not symmetrical about @ = @r. For the special cas | 
infinite at the leading edge, 9 7 + - 6;, as in 8 = 0, however, this symmetry is present and the} 
the case of the isolated airfoil. The integration covers analogs of Eqs. (3.21) and (3.22) can be obtained 
From Eqs. (4.25) and (4.26), it follows, by separating | 
odd and even parts, that, for camber, 
be 
7,/(0) = (1 [G(¢) G(@)] [sin (cos — cos dod 
0 
} Pe 
u,'(0) = —(1, IG(o) | [sin ¢/(cos — cos 6)] + [K*? — cos? [sin ¢/(1 — cos ¢)]} do 
0 
7 se 
= —(4K*a K) G()2,'(@) cot (6/2) dd = ze 
0 0 of 
while, for thickness, [, = 0 and 4) 
e\ 
v,/(0) = an) f IG(@) G(@)| [sin ¢/(cos — cos 4)| dd 
0 
u,'(@) = —(1 IG(@)/G(@)] [sin 6/(cos — cos dg 
0 (I 
Here x = (r/m) tanh! (cos 0/K), & = (r/m) tanh! (cos ¢/K) 
while dé = —2aG(¢) do 
al 
with G(¢) = (K 7/27a) sin ¢/(K? — cos? ¢) as 


i 

3 
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le know The equations for camber become 
v./(x) = (1/7) | coth r) — x)dx (4.51a) 
| 
conditio 11 — K tanh (ax/r) 1 1 + A tanh (7é/7) l 
= \ - coth — (¢ — x) — (4.31b) 
rted int; 1+ A tanh 7 2 1 — K tanh (2é/7r) T K 
pressions 
=2 u,(&) dé (4.31¢) 
) (4.99 The thickness equations are 
= — = coth — — x) dé (4.32a) 
) T J —c/2 1 — tanh 7) T 
u, (x) = v,/(~) coth — (£ — x) dé (4.52b) 
T J -—c/2 
(4.30a 
Eqs. (4.3la) and (4.32b) are Schlichting’s® Eqs. (5) and (6) for zero stagger with = 2u,’(&), g(g) = 22," 
It is worth noting the modifications of Eqs. (4.30) and (4.31) employing the classical boundary condition. For 
camber, 
(4.30b = I cos a(dy, dx) — V sina 
so that 
(4.30¢ 1 — A tanh (wv 7) jA* l'cosa dy, 1 + A tanh 7) 
j 1+ A tanh (ax 7) UK T J » dé 1 — K tanh (7é 7) 
n (07 + coth — — «) — K dé (4.53a) 
T 
Vv shape 
nd cam: | sin a K* V cosa 2 dy, 1 + K tanh 7) (4.331 
K K J Ge 1 — AK tanh (rt 7) 
also, 
funda- M./pV = (2Vr/x) sin a In (AK /A*) — cosa (dy, dé) cos~' [(K* K) cosh (ré 7) (4.33¢) 
| —c/2 
ted air- | 
(43 For thickness, 
ial case | c/2 Klingeman in Eq. (64) of reference 10. As far as is 
nd the > cos @ 7) (dy, dé) coth known, the expression for — upp is new. 
tained. | — x)dé (4.33d) 
‘ (5) NUMERICAL SOLUTION OF THE INTEGRAL EQUATION 
Eqs. (4.33a), (4.33b), and (4.33c) have essentially 
been given before by Dorr* and Nickel." Eq. (4.33d), While it is possible to solve explicitly the direct prob- 
} appears to be new. lem for the airfoil in cascade using the boundary 


Parabolic Arc Airfoil 


do 
As an example of the use of the formulas of the last Z- PLANE 


section, application is made to the parabolic airfoil of 
zero thickness shown in Fig. 4.4. If h is the height 
of the are at mid-chord, the equation is y = (1 — 


The integrals appearing in (4.29) can be 


evaluated. Here 
ut — urp) (V cos a) = (Shr we?) X 
lcos 8 are tan (sin @ A*) — sin 8 are tanh (cos 6 K)] 


— Trp)/(V cos a) = 

(16h (7 c)? In (K A*), — = O 
These results agree with the known ones for the same ee Gm, 
airfoilask—~>Owith 8 = 0. The circulation is the same -Cc/2 0 c/2 
as that found for a similar velocity distribution by Fic. 4.4. Parabolic are cascade. 


| 
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pressure 
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condition 


v= u(dy dx) 


it is considerably easier to carry through the solution numerically. 


OF THIN 
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on the chord 


The integral equation (4.28) 


1 (@ — 6r) u(@) dy 
= — | J cot (A? — cos? 6) — cot (K? — cos? 6r) + 
Jy (6 — Or) 2 fo dx 
jAA* sin a + AJ sin cos (a + 8) — A*J cos (a + 
J? sin (0 — Or) 
is solved by the usual iteration procedure. A given 


approximation to u(@) is inserted in the integral on the 
right and yields the next approximation u(@). If the 
first approximation is taken to be u(@) = I’ cos a, the 
second approximation will be the solution associated 
with the classical thin airfoil theory. Two or three 
further approximations will usually suffice for con- 
vergence. The solution so obtained is referred to as 
the “integral equation solution.”’ 

For the reader who wishes to solve Eq. (4.28) using a 
desk calculator, it is suggested that 10 or 12 points be 
used at equally spaced values of 6 (67 and 67 + 7 are 
avoided because of the complexities at these poitts). 
This choice will give a natural condensation of points 
on the airfoil near the leading and trailing edges which 
is desirable. The equal spacing has the additional 
advantage of simplifying significantly the trapezoidal 
For 2.V equally 


approximation to the integral.’ 


spaced points on the interval 0, 27, 


(t 27) cot } [(mm N) — 2| do = 
J 


I<P<N-1 


ph odd 


(1 NV) cot (px 2N) f[m + p N)] - 


flm — p(w N)]f 


The remaining expressions within braces are also com- 
puted independently of the airfoil contour under study. 
Note further that it is necessary to solve the integral 
equation for only two values of a in order to obtain the 
complete dependence upon a@ [ef. Eq. (3.27)]. It is to 
be remarked finally that average properties of the ve- 
locity distribution such as turning angle can be deter- 
mined fairly well with a single approximation. 

These numerical solutions have been obtained at the 
United Aircraft Corporation Research Division for ex- 
tensive ranges of each of the parameters involved. 
The abscissa x instead of 6 was used as the independent 
variable in order to obtain pressures at prescribed 
locations without interpolation. The interval of inte- 
gration was divided into 54 intervals, the gage of the 
subdivision being finer at the leading and trailing edges 
of the airfoil. The work was performed on an IBM 
Card-Programmed Calculator and later on a 704 
Electronic Digital Computer. In its final form on the 
‘04 the entire calculation—which includes iteration 
until a satisfactory degree of convergence is obtained 
lor the velocity distribution, determination of pressure 
distribution, turning angle, inlet and exit angles, and lift 
coeflicients for a wide range of operating conditions 


requires approximately one minute. 


(§) COMPARISON WITH Exact SOLUTION 


For the purpose of indicating the degree of accuracy 
of the integral equation solution developed here, a 
comparison has been made between our calculated 
pressure distributions and some exact theoretical data 
published in reference 8. Results of this correlation 
are shown for the case of an isolated airfoil in Fig. 
6.1 and for the case of a cascade of airfoils in Fig. 6.2. 
Theoretical distributions calculated using 
the boundary conditions of classical thin airfoil theory 
(referred to as “classical solution’’) have been super- 
posed on these plots also for purpose of comparison. 
As illustrated, the integral equation solution gives 
pressure distributions which agree quite satisfactorily 
with the exact results. Of particular significance is the 
finite pressure peak at the leading edge in contrast to 
the characteristic infinite value for the classical thin 


pressure 


airfoil solution. 


(7) COMPARISON EXPERIMENT 


A comparison of experimental and calculated data 
has been made to demonstrate the adequacy of the 
present theory in simulating the physical model. Typi 
cal results presented in Fig. 7.1 for a staggered cascade 
of low cambered NACA 65(4)10 airfoils show the 
agreement between theory and experiment’ to be 
excellent. It is of interest to note that with the ex- 
ception of the leading-edge region of the airfoil, the 
classical pressure distribution superposed on Fig. 7.1 
is also in close agreement with experimental data. 

Although the nature of approximation employed in 
establishing Eq. (4.28) makes use of the assumption of 
thin, slightly cambered airfoils, extensive experience 
has indicated that reliable pressure distribution infor- 
mation is obtained even for highly cambered profiles 
of the type presently used in axial flow turbomachines. 
As examples, Figs. 7.2 and 7.3 are included to illustrate 
the good agreement obtained between theory and ex- 
periment for cascades of NACA 65(12)10 and 65(18)10 
series airfoils, respectively. 

While knowledge of pressure distributions is poten- 
tially useful to evaluate loss coefficients, boundary- 
laver development, and flow separation characteristics 
of airfoils in cascades, a comprehensive knowledge of 
flow deflection angles is also required to provide proper 
matching of compressor stages. In Fig. 7.4, a com- 
parison between calculated and experimentally meas 
ured turning angles is presented which gives some indi- 


| | 
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quence allows for substantial reductions in experj.| 


NACA 65(18)10 AIRFOIL mental programs required to supply these data. a 


B =68° To demonstrate the flexibility and usefulness of the = ™ 
CH present theory for this purpose, calculations wer " 
-20 performed to obtain the cascade performance char. © 
INTEGRAL acteristics of NACA 65 series airfoils. The NACA @ 
- SOLUTIO series section was selected for this study for several . 
S -15 |, CLASSICAL SOLUTION -- reasons. 
The 65-series section is representative of a class oj : 
i or DATA ° airfoils which have good high Mach Number perform. t 
-1.0 ance characteristics and which, as a consequence, ar 
oO being designed into present and advanced turbojets a 
iL ~~ The availability of extensive experimental data on this) 4 
Le series of airfoils permits comparison and _ correlation 
! -05 between theory and experiment. The simple empirical ¢ 
O rule,'* for flow turning-angle, referred to as Carter's 
< o\ or Howell's rule, cannot be directly applied to the 65 I 
3 e) Nc series type of airfoil because of the arbitrary manner P 
= oN\ of defining the camber angle for mean camber lines } ‘ 
n which have infinite slopes at the leading and trailing 
aw ? —— An illustration of the effect of variation of blade a 
a camber angle, thickness-chord ratio, and gap-chord d 
ratio on the flow turning characteristics of the NACA 
a 25 50 75 100 65 series airfoils is shown in Figs. 8.1, 8.2, and 83.) _ 
PERCENT CHORD, x/c 
Fic. 7.3. Comparison of theoretical and experimental pressure 
distributions on airfoils in cascade. NACA 65(4)I0 AIRFOIL 
Bou = 56.5° 
cation of the accuracy of the integral equation solution t/c=1.0 ; 
in regard to predicting air turning angles for cascades. 20° . : 
The good agreement between theory and experiment EXPERIMENTAL ‘ 
obtained in this specific case is typical of that obtained DATA o 
when the correlation is made in the unstalled operating i 
range of the cascade configuration. THEORY —— 
15° 
(8) APPLICATION OF INTEGRAL EQUATION SOLUTION 
IN COMPRESSOR DESIGN (3) 
In the initial phase of compressor design, the flow ws 10° | 
velocities and angles at the inlet and exit of the several ) ¥ 
blade passages are usually specified from prior con- = ) 
siderations, such as allowable pressure rise per stage, 
critical Mach Number, and fluid losses. The blade © 
designer must then select a cascade, from a multitude = 5° 
of possible combinations, to satisfy these initial flow a 
requirements. To accomplish this objective, com- = 
prehensive information must be provided which ade- 
quately define the effect of variation of the several 0° 
cascade parameters—e.g., gap/chord, stagger angle, ; 
blade profile—on cascade flow characteristics. In the 
absence of a convenient and reliable theoretical analysis, 
these data are obtained by extensive experimental tests | 
and the development of semiempirical correlation -5° } 
schemes.'*! The simplicity of application of the pres- 30° 40° 50° 60° 
ent integral equation solution makes it possible to INLET ANGLE, B, r , 


supplement rather rapidly available information on 
PI ada F Fic. 7.4. Comparison of theoretical and experimental flow turn- 


cascade performance characteristics, and as a conse- ing angle. 
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respectively. The results presented in these Figures 
are typical of those obtained in our studies over a wide 
range of cascade geometries. The fact that the turn- 
ing angle is very nearly a linear function of blade camber 
9*, gap-chord ratio 7 ¢, thickness ratio ¢ c, stagger 
angle Ben(= 90° — @), and inlet angle 3; over the 
range 0 < O* < 40°,05 15,0 < 0.12, 
10° < Bey S 60°, 40° < 2B; < 60°, leads to the formu- 
lation of obvious empirical rules for the prediction of 
the cascade performance of this airfoil series. 

In conelusion, it should be noted that if consider- 
ation is limited to thin airfoils having small camber, 
and if definition of the pressure peak is mot required, 
then the classical cascade solution gives flow turning 
characteristics that are in very good agreement with 
results obtained using the integral equation solution. 
This observation is based on a fairly extensive com- 
parison of theoretical data for cascade geometries cover- 
ing the range of 0.53 < 7 ¢ < 1.50 and 30° < Bey < 60°. 
The maximum difference of flow turning angles between 
the two solutions exists at the low gap-chord condition, 
and, as illustrated in Fig. 8.4, is generally less than 1 
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APPENDIX 


Proof of Fundamental Principle and Its Corollary 


It may quite properly be said that the roles of the 
Fundamental Principle and its Corollary might have 
been interchanged, for it is convenient to prove the 
Corollary first and to obtain the other as a consequence 

Assume that the function /(¢) is analytic in ¢ <a, 
Then by the Cauchy Integral 


continuous on ¢ = a. 
Formula, for = <a, 


= (C1 2212) [MOde — 
=a 


In particular, 


fO) = ani) de 


By the Cauchy theorem, since /(¢) [€ — (a* ¢)] is an- 


alytic inside ¢| = a,so that 


(1 ani) de — (a? = 0 
=a 


it follows that 


degree. 
NACA 65(XX)OO AIRFOILS NACA 65(2.8)09 AIRFOIL 
GAP/ CHORD = 1.070 
60° Bou = 50 
30°; ' 
® | 
oO 
8° 
© B,= 45° 
4 
B,= 60° — 
a 
4 = B,= 50° 
| O 04 08 1.2 1.6 
20° 30° 50° GAP/ CHORD, 
a Fic. 8.3. Effect of gap/chord variation on flow turning angle 
BLADE CAMBER ANGLE, @ 
Fic. 8.1. Effect of variation of blade camber angle on flow turn- 
ing angle. 
NACA 652.8(XX) AIRFOILS 
GAP/ CHORD = 1.070 NACA 65(2.8)09 AIRFOIL 
= 60° GAP/ CHORD = 0.536 
30° T 24°, 
| INTEGRAL EQUATION Bcou= 60 
SOLUTION —— 250° 
° @ CLASSICAL SOLUTION —— IN 
|B, 40 + Bou 40 
Bou * 30° 
2 Zz ° 
> 
B,= 60° 
fe) 2 4 6 8 10 12 o° 10° 20° 30° 40° 50° 60° 
BLADE THICKNESS, t/c INLET ANGLE, B, 
FiG. 8.2. Effeet of variation of blade thickness on flow turning Fic. 8.4. Comparison of turning angle using classical and in- 
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tegral equation solution. 


j 
| 
| 
5 
16) ‘ 
— 


9 


= (1 2m) | f(ae’’) [aX(e""? — 
0 


The limit may be taken as r > a 


fiae") = f(0) + | de 


= f(0) — f(ae’®) cot do 
o 2 


This singular integral is to be interpreted from the 
point of view of Cauchy Principal Values: for 0 
<0 < Se, 


cot — do = 
0 


2r 
lim + ) fae) cot dg 
0 


The case 6 = 0 or 6 = 2m is handled by a shift of 
interval of integration. 
The Fundamental Principle itself may be obtained 


from an application of the transformation Z = a® ¢. 
Let f(¢) = F(Z) = Then 


2n 
= F() — | F(ae~'®) cot 
0 


Replacing @ by —6@, and ¢ by 27 — @, and using the 
periodic property of F as a function of 6, gives 


F(ae") = F(~) +. F(ae’®) cot 
0 2 


do 


REFERENCES 


' Lighthill, M. J., 4 New Approach to Thin Airfoil Theory, 
Aero Quarterly, Vol. 3, pp. 193-210, November, 1951. 

2Van Dyke, M., Subsonic Edges in Thin-Wing and Slender- 
Body Theory, NACA TN 3348, 1954. See also, Van Dyke, M., 
Second-Order Subsonic Airfoil-Section Theory and Its Practical 
Application, NACA TN 3390, 1955. 

3 Pistolesi, E., On the Calculation of Flow Past an Infinite 
Screen of Thin Airfoils, NACA TM 968, 1941. 

4 Lieblein, V., Zur Berechnung der Auftriebscharacteristik eines 
Profils im Gitterverband, Ingenieur-Archiv, Vol. 18, No. 5, pp. 
281-290, 1950. 

5 Scholz, N. On the Extension of Glauert’s Theory of Thin Air- 


JOURNAL OF THE AERO/S 


ti=a 


PACE SCIENCES—JULY, 1958 i 


= (1/27) f(ae’®) \ [1 (ae’® — + [1 (ae’® — r)] — (1 /ae’*)! ae’ do 
0 


210 o » 16 
(ae’® — re") (ae’® — ate” r)| 


foils to Profiles in Cascade, Journal of the Aeronautical Sciences 
Vol. 18, No. 9, pp. 637-639, September, 1951. See also Sty; 
mungsuntersuchungen an Schaufelgittern, V.D.1. Forschungsheft 
442, 1954. 


6 Schlichting, H., Berechnung der reibungslosen inkompressibley dy 
Strémung ftir ein vorgegebenes ebenes Schaufelgitter, VDA. For. | ye 
schungsheft 447, 1955. See also Proceedings of the Conference re 
on High-Speed Aeronautics, Brooklyn Poly. Inst., New York, sat 
1955. bes 

7 Tyler, R. A., The Available Theoretical Analyses of 
Dimensional Cascade Flow, Nat. Res. Council of Canada, Aero as 
Note AN-4, 1949. | 

8 Garrick, I. E., On the Plane Potential Flow Past a Lattice of | ne 


Airfoils, NACA TR 788, 1944. : 
Dorr, J., Strenge Lésung der Integralgleichung fiir die Stré- 

. . . 

mung durch ein senkrechtes Fligelgitter, Ingenieur-Archiv, Vol. 19, 


pp. 66-68, 1951. 

” Nickel, K., Zusats su J. Dérr: Strenge Lésung der Integral- 


gleichung fiir ein Fligelgitter, Ingenieur-Archiv, Vol. 20, pp 
6-7, 1952. See also Uber spesielle Systeme von Trag fliigelgittern Be 
II, Ingenieur-Archiv, Vol. 23, pp. 107-108, 1955. 

1 Durand, W. F., (Ed.), ‘‘Aerodynamic Theory,’”? von Kar- 
man, Th., and Burgers, J. M., General Aerodynamic Theory 
Perfect Fluids, Vol. I1, Julius Springer, Berlin, 1934. Also, 
Durand Reprinting Committee, Pasadena, California, 1943. 

2 Allen, H. Julian, General Theory of Airfoils Having Arbitrar) 
Shape or Pressure Distribution, NACA TR 8833, 1945. 

18 S6hngen, H., Die Lésung der Integralgleichung und deren 
Anwendung in der Trafliigeltheorie, Math Zeitschift, Vol. 45, pp 
245-264, 1939. 

4 Klingemann, G., Verfahren zur Berechnung der theoretischen 
Kennlinien von Turbomachinen, Ingenieur-Archiv, Vol. 11, pp. 
151-178, 1940. 

% Smith, C. B., A Solution for the Lift and Drag of Airfoils with 
Atr Inlets and Suction Slots, Journal of the Aeronautical Sciences, 
Vol. 16, No. 10, pp. 581-589, October, 1949. ' 

Watson, E. J., Formulae for the Computation of the Functions 
Employed in Calculating the Velocity Distribution about a Given 
Aerofoil, British R&M 2176, May, 1945. 

7 Herrig, L. J., Emery, J. C., and Erwin, J. R., Systematic Two- 
Dimensional Cascade Test of NACA 65-series Compressor Blades 
at Low Speeds, NACA RM L51G31, September, 1951. 

8 Howell, A. R., Fluid Dynamics of Axial Compressors. De- 
velopment of the British Gas Turbine Jet Unit, Inst. Mech. Eng. | 
London, England, ASME report, 1947. 

18 Sorensen, H. A., Gas Turbines, Ronald Press Company, New 
York, 1951. 


ology, etc. 


1958. 


IAS Twenty-Seventh Annual Meeting r 
Hotel Astor, New York City, January 26-29, 1959 


The scope of the meeting will be broad to cover the many fields of interest to IAS members. The meetings com- a 
mittee plans to schedule papers on aerodynamics, hypersonics, aeroelasticity, missiles and rockets, missile propul- 
sion, space propulsion, satellites and spacecraft, space medicine, VTOL/STOL aircraft, rotary wing aircraft, meteor- 


Members or organizations wishing to have papers considered for presentation at this meeting should submit out- \ 
lines or short abstracts to the Meetings Committee, 2 East 64th St., New York 21, N.Y., no later than August 1, : 


424 
ve 
2G 
£ 


il Sciences 


also Sty 


chungsheft 


Npressible n 


For. 


Conference 


ew York 


Two- 


ida, Aero 


s of 


Lattice of 


die Stré- 
Vol. 19, 


I ntegral- 
20, pp 


'gelgittern 


on Kar 
Theory 
1. Also, 
143, 

1 rbitrary 


deren 
45, pp 


retischen 
11, pp 


with 


ciences, 


inctions 
Given 


ic Two- 


Blades 


. De 
Eng 


New 


A Rational Determination of Loads and Exit 
Velocities of Cartridge Ejected Stores 


H. WOLF* anp S. PINES** 


Republic Aviation Corporation 


SUMMARY 


A rational analytical procedure is presented which enables the 
jynamics engineer to determine peak gun forces, ejection time 
histories, separation velocities, and dynamic structural amplifi- 
cation factors of cartridge ejection mechanisms for store release. 
Specific computations are carried out illustrating the effect of 
various pertinent parameters such as initial chamber volume, 
power piston cross sectional area and powder burning rate on the 
store exit velocity, peak pressure, and gun loads 

The paper derives the thermodynamic equations of the chemi- 
cal energy exchange in a form convenient for solution of the 
problem. Since the pressure time history of the burning powder 
isa sensitive function of the chamber volume, the motion of both 
the ejected store and the breech reaction are required. The 
equations of motion of the elastic restraining structure, and the 
ejected store, are derived using standard procedures. The re- 
sulting system of nonlinear simultaneous differentia! equations 
is presented. 

The resulting solutions, obtained by analog computing pro- 
cedures, are presented in the form of curves which may be used in 


ejection gun design. 


SYMBOLS 


gas pressure 
= gas volume 


7 = absolute temperature 
E = internal energy of gas 
Q = heat introduced by combustion 
K = heat equivalent of burning powder mass 
A = piston area 
\ = relative vertical displacement of piston to gun 
ds = vertical store displacement 
dg = vertical gun displacement 
F = peak gun reaction force 
w = natural frequency of gun support structure 
Ms = store mass 
m; = generalized mass of gun and supporting structure 
d/dt = derivative with respect to time 
i = initial chamber volume 
= adiabatic constant of combustion gases 
( = burning rate constant 
R = gas constant (per unit mass) 
m = mass of gas 
c = specific heat at constant pressure 
t. = specific heat at constant volume 
p = density of powder 


INTRODUCTION 


_ CHOICE of a gun for store ejection is based on 
many practical details. In modern aircraft, where 
ejection is often required at high Mach Numbers and 
at high indicated speeds, the problem of separation and 

Presented at the Aeroelasticity—I Session, Twenty-Sixth 
Annual Meeting, IAS, New York, January 27-30, 1958. 

* Group Computing Engineer. 

** Chief Flutter and Computing Engineer. 


high dynamic wing reaction loads makes this choice 
particularly critical. Consequently it is considered 
timely to give gun designers and dynamics engineers a 
rational means of predicting exit velocities and peak 
loads as a function of some of the pertinent parameters. 

Many dynamics engineers attempt to solve this prob- 
lem assuming the time history of the explosion as 
known, going to test results of rigidly mounted gun 
discharges for typical time histories. Unfortunately, 
this procedure is doomed to failure since the pressure 
time history is a sensitive function of the instantaneous 
chamber volume which in turn is a function of the mo- 
tion of both the store and the wing gun structure. Thus, 
the pressure time history of a rigidly mounted gun is 
different from that of a flexibly mounted gun. 

The equations of motion of the store and wing gun 
are well known, and by adding the thermodynamic 
pressure equations derived herein, a rational unified 
attack of the problem becomes feasible. 


THE THERMODYNAMICS OF THE GUN 


In its simplest form, the store ejection gun consists of 
a chamber closed at one end by a movable piston to 
which the store is attached (see Fig. 1). A propellant 
is burned in the chamber, generating a gas and thermal 
energy. The pressure of this gas on the piston produces 
the ejection force, and the reaction to this force pro- 
duces the load on the gun supporting structure. In 
addition, the gas pressure determines the burning rate 


of the propellant. Although the pressures within the 


qs 


Ms 


Schematic of gun store and flexible wing. 
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gun are usually quite high, the ideal gas laws hold to a 
sufficiently high degree of approximation, so that the 
following theory, which is based on these laws, permits 
predictions correct within the experimental error. 
In order to account for the variable mass of the gas 
present, the ideal gas laws were used in the form® 


pv = RmT (1) 
dE dt = (dQ dt) — p(dv dt) (2) 
E = (3) 


If dQ dt is taken as zero, the above three equations re- 
duce to the adiabatic gas law. However, in our case, 
the burning powder introduces heat into the system. 
The amount of heat introduced is proportional to 
the mass of powder burnt, if heat losses are neglected. 
The time involved is usually so short that this may 
This results in the equation 


dQ dt = K(dm dt) (4) 


safely be done. 


The mass rate of burning is primarily a function of the 
chamber pressure (for a given cartridge), so that 


dm, dt = f(p) (5) 


Various forms have been suggested for the function /(/) 
and have been tried. DeSaint Robert’s law 


f(p) = ep" (6) 
with » = 0.7 approximately, has been found most satis- 
factory.! 

Eliminating /, m, Q, and 7 from Eqs. (1) through (6) 
and use of the relations 


C,/C, 


results in 
v(dp dt) + yp(dv dt) = cK(y — 1)p*" (7) 
After the powder is exhausted, Eq. (7) becomes 


(7a) 


v(dp dt) + yp(dv/dt) = 0 


which is the differential form of the adiabatic law. 
Integration of Eq. (7) gives the pressure as a function 
of time if the time history of the volume available to 
the gas is known. 

Before proceeding to the determination of the equa- 
tions required to describe the time history of the cham- 
ber volume, it is convenient to indicate the procedure 
for evaluating the constant cK (y — 1). The method 
proposed herein is to use the initial pressure time history 
of a typical exploding charge. Initially, and for a 
short time interval following, the rate of change of 
volume dv dt can be taken to be zero and the equation 
may be integrated independent of the store and wing- 
gun motion. The solution is 


— = 0.3cK(y — (8) 
Using this relation, an accurate determination of 
cK (y — 1) is possible directly from test data of a given 
charge. Typical values are entered below for several 
types of cartridge combinations: 
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to 


Mk. II 0 
Mk. II 0-1” Navy and 1 OM 2” 
cK (y — 1) = 7460 
Mk. II O— 1” and 1 Mk. 1 
cK (y — 1) = 7930 


_ 


Cartridg, 


I] 3” Navy Cartridg 


THE EQUATIONS OF MOTION 


The rate of change of volume consists of two parts 
(a) the volume increase due to the displacement of the 
piston relative to the gun, and (b) the volume vacated 
by the burning powder. 

The second part is small compared with the first after 
the explosion has been in progress for some time and the 
piston has attained an appreciable velocity. Its omis- 
sion, however, would cause the pressure to rise too rap 
idly in the initial stages. 
therefore, 


The volume equation js 


dv dt = A(dx dt) + (c p)p*? 4 


Eq. (9) is given in a form that remains valid even if the 
piston cross section changes, as in a telescoping gun. 

dx dt in Eq. (9) is the relative velocity of gun and 
store and depends on the motion of both; thus, 


dx dt = (dqs dt) — (dq« dt) (10 


The equations of motion for the gun and the store 
have to be used to obtain gs and gg. For the simple 


model (see Fig. 1) these equations are 
ms(d*qs dt?) = Ap (11 


me(d°*qq dt?) + mgwqg = —Ap (12 


In the special case of a rigidly supported gun, Eq. (12 
is replaced by 


Gg = 


The freely recoiling gun can be treated as a rigid gun 
with equivalent store mass given by msmg (ms + me 
the harmonic mean between gun and store mass. 

In the practical case of an ejection gun mounted on 
an aircraft wing, Eq. (12) has to be replaced by the 
Lagrangian equations of motion. In this analysis the 
normal modes of the wing structure and the generalized 
masses associated with them are used. The contribu- 
tion of each mode to the relative velocity dx dt has to 
be found. In general, dx ‘dt will become a linear com- 
bination of all the modal velocities. The pressure as 
well as aerodynamic and friction forces are introduced 
by the work done in each mode. A detailed discussion 
of the derivation of these equations, especially the 
manner in which aerodynamic and friction forces are 
introduced, is beyond the scope of this paper and will 
be the subject of a future note. 


DISCUSSION 


Eqs. (7), (9), (10), (11), and (12) were integrated 
simultaneously on an analog computer. Fig. 2 com- 
pares pressure-time curves obtained by Chicago Mid- 

yay Laboratories? with curves computed on the basis 
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of the above theory for a number of stores and initial 
chamber volumes. The agreement is seen to be good. 

Exit velocity and peak load on the gun supporting 
structure are the most important results to be obtained. 
These quantities have been plotted for a 2,000-lb. store 
as a function of the piston area, the initial chamber 
volume, and the mass burning rate. It is desirable to 
achieve high exit velocities without the use of exces- 
sively high peak forces. As might be expected, a higher 
peak gun load will in general result in a higher sepa- 
ration velocity. Fig. 3 shows the dependence of exit 
velocity and peak force on the initial volume and the 
mass burning rate for constant piston area. In order 
to obtain a higher exit velocity, the mass burning rate 
must be increased or the initial volume decreased. 
From the relative slope of these curves it is seen that 
increasing the mass burning rate will give a gain of exit 
velocity at lower cost in peak force than a decrease in 
initial volume. Figs. 4 and 5 show the dependence of 
exit velocity and peak load on the cross-sectional area 
of the power piston. It can again be seen that decreas- 
ing the area of the piston will increase exit velocity at 
smaller cost in peak force than either a decrease in initial 
volume or an increase in the burning rate. This is due 
to the fact that a small cross-sectional area makes the 
rate of change of volume dv dt small, even for high store 
velocities, and thus prevents a rapid drop-off of pressure 
after the peak has been attained. The required strength 
of the piston, however, places a lower bound on its 
cross-sectional area. 

Fig. 6 shows the effect of flexibility of the gun sup- 
porting structure on loads and velocities, a typical wing 
mass being used. It shows that a more flexible support 
will materially lower the loads on the gun and also the 
absolute store velocity. The relative velocity of store 
and gun, however, will increase if the stiffness and nat- 
ural frequency are lowered below a critical value. Past 
this point, separation from the wing will actually be 
improved by increasing flexibility of the gun support. 
It should be noted, however, that if clearance from the 
aircraft, rather than just clearance from the wing, has 
to be considered, the store velocity dgs dt rather than 
the separation velocity dx dt is the important param- 
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eter, and, for very low wing rigidities, the gun velocity 
accounts for the major part of the separation velocity 
It is also quite apparent that loads and velocities are 
quite different for the flexibly supported gun, so that 
an analysis using pressure time curves obtained from g 
rigidly mounted gun will furnish incorrect results. 

Lengthening the power stroke of the gun furnishes g 
method of increasing exit velocity at no cost in guy 
loads, if the friction of the piston in the barrel can be 
controlled. Long pistons, however, will make the 
moment due to the air drag forces large, with an attend- 
ant rise in friction. Consequently, little exit velocity 
is gained during the latter part of the stroke. 


CONCLUSIONS 


The method presented adequately predicts gun loads 
and exit velocities. 

Flexibility of the gun supporting structure materially 
alters loads and exit velocities. ; 

Decreasing the piston cross-sectional area increases 
exit velocities with small increase in gun loads. This 
area, therefore, should be taken as small as _ possible, 
consistent with the required piston strength. 

Ejection velocities may be increased by decreasing 
the initial volume of the chamber. This is accompanied 
by a notable increase in peak gun loads. 
velocities also may be increased by using a cartridge 
with higher mass burning rate, accompanied again by 
increasing peak loads. 

Wing flexibility reduces gun loads and store velocities. 
It will, however, increase separation velocities if it is 
large enough to permit considerable wing motion. 


Ejection 
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A Critical Strain Approach to Creep 
Buckling of Plates and Shells’ 


GEORGE GERARD* ann ARTHUR C. GILBERT** 
New York University 


SUMMARY 


A creep buckling hypothesis presented previously suggested 
that the known solutions for short-time buckling can be used in 
terms of critical strain to estimate creep buckling. 

The major portion of the paper is concerned with a creep buck- 
ling investigation of aluminum cylinders under compressive and 
The test techniques and equipment are described 

The results of the experiments are correlated 


torsional loads 
in some detail. 
with the critical strain approach, and good agreement is obtained 
both in terms of critical strain and critical time. Some _ prelimi- 
nary data on creep buckling of cylinders under cyclic load condi- 
tions are also presented. 

An examination of test data on creep of long flat plates under 
compressive loads indicates that the theory may be adequate for 
plates of small width-thickness ratio. Plates of large width- 
thickness ratio which can buckle upon the initial application of 
load require an analysis of the post-buckled state under creep 


conditions. 


SYMBOLS 


h = plate width, in. 

E = modulus of elasticity, psi 
E, = secant modulus, psi 

G = shear modulus, psi 

h = plate thickness, in. 

L = length, in. 

R = radius, in. 

t = time, hours 


= shear strain, in. /in. 


€ = strain, in./in. 

vy = Poisson’s ratio 

o, = applied creep stress, psi 

tT, = applied creep shear stress, psi 
T-- = critical shear stress, psi 


(1) INTRODUCTION 


[ REFERENCE 1, an hypothesis was advanced for es- 
timating the critical time for creep buckling of col- 
umns, plates, and shells. This hypothesis suggested 
that, as an engineering approximation, creep buckling 
occurred when the critical strain determined from elas- 
tic and plastic stability considerations was reached in 
the creep process. 
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The critical strain approach indicated reasonable cor- 
relation with some creep buckling test data on columns 
as indicated in reference 1. In order to obtain basic 
creep buckling data on structural elements other than 
columns, an investigation was conducted on cylinders 
under compression and torsion. Because most of the 
equipment was specially designed for this investigation, 
the principal features and techniques are described be- 
low. 


(2) EXPERIMENTAL EQUIPMENT AND PROCEDURES 


Material 

The material selected for this investigation was com- 
mercially pure 3003-0 aluminum which is not subject 
to the structural changes associated with the heat-treat- 
able alloys at elevated temperatures. To minimize 
variability in mechanical properties, all specimens used 
in the investigation were machined from a single 16 ft. 
length of extruded tubing. The tubing was circular 
in shape with an inner diameter of 1.045 in. and 0.135 
in. wall thickness. 


Creep Units 

The compressive creep unit with its associated equip- 
ment is illustrated in Fig. 1. It consists of a conven- 
tional loading frame with a pneumatic load cell mounted 
in the base. The load cell contains an essentially fric- 
tionless seal, and, therefore, the regulated nitrogen tank 
pressure supply used for loading is monitored at the 
cell for load indication within +1 2 per cent accuracy 
With this system, a constant load can be maintained for 
an indefinite period. 

The torsional creep unit shown in Fig. 2 utilizes a dead 
weight attached to a wheel for torque application. 
The movable loading head fits into a carefully machined 
square hole in the wheel which in turn is supported 
on roller bearings attached to the frame of the machine. 
A counterbalanced shaft, centered through the wheel, is 
pinned so as to permit axial motion of the specimen. 
The counterbalance feature prevents bending of the 
specimen. Calibration of the machine indicated that 
friction in the system was less than 0.8 per cent of the 
minimum load applied in the creep program. 

At the fixed end, the applied torque is reacted by a 
hollow shaft firmly attached to the frame. The fur- 
nace control or sensing thermocouple is passed through 
this bar and is attached directly to the specimen. 

Both ends of the torsion specimens have threaded 


4 
| 
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Fic. 1. Compressive creep unit. 


portions sufficiently long so as to permit the threads to 
bottom firmly in the grip. At elevated temperatures, 
differential expansion between the steel grips and the 
aluminum specimen firmly locks the specimen in the 
grips. 


Extensometers and Recorders 


In the compression tests, specimens are positioned in 
the furnace between the upper and lower plattens. 
Creep or axial deformation is measured by observing 
the motion of the loading platten with respect to the 
stationary platten, as multiplied through the extensome- 
ter linkage shown in Fig. 1. 

For short time tests, the deformation is magnified by 
the extensometer linkage and read on a dial gage capable 
of being read to 0.0005 in. With a 10:1 magnification 
of the extensometer, the deformation can be measured 
within 50 microinches. 

Long-time creep deformation was measured by re- 
cording the motion of the loading platten with respect 
to the stationary platten, as magnified through the ex- 
tensometer linkage. Wire rope of approximately 0.018 
in. diameter high strength steel is anchored to the ex- 
tensometer sector and passed through a well lubricated 
casing to a Foxboro temperature recorder, reworked to 
record motion rather than temperature. 

By adjusting the lever system arrangement within the 
recorder it was possible to magnify the motion and 
maintain linearity over the range of recorder operation. 
The recorder-extensometer system was calibrated 
against the motion of the loading platten by use of an 
highly accurate optical system. It was determined 
that the recorder charts could be read to within an 
accuracy of 400 microinches. The clock mechanism 
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and the lever system of the Foxboro apparatus pro- 
vided the required timing and writing mechanism, 
Maximum error in the timing device was 0.5 per cent. 

For short-time torsion tests, a dial gage capable of 
being estimated to 0.0002 in. was positioned with re- 
spect to the loading wheel such that the wheel rotation 
and consequently shearing strain was measured. The 
torsional creep of the specimen was also obtained from 
the rotation of the loading wheel. The rotation was 
transmitted through the supporting rollers to a pulley 
mounted concentrically with one of the rollers. Wire 
rope transmitted the motion of the pulley to a re- 
worked Foxboro recorder utilizing an arrangement sim- 
ilar to that used in the compression creep machine. 

The recorder was calibrated against the motion of the 
loading pulley by affixing an accurate protractor to the 
loading pulley. With a filar eyepiece focused on the 
protractor, it was possible to measure pulley motion to 
0.01°. The resulting calibration error in the recording 
system was of the order of 0.2° or a shear strain of 
approximately 500 microinches. 


Heating Equipment 


Both the compression and torsion creep units were de- 
signed to use a Marshall tubular furnace having nine 
heating zones. Consistent creep results require a con- 
stant temperature along the specimen gage length. 
Therefore, a specimen with thermocouples located at 
three points along the gage length was installed in the 
machine simulating a creep test and the current to each 
of the heating coils was then adjusted so that the max- 
imum temperature variation along the gage length was 
2'/,°F. 


Fic. 2. Torsional creep unit. 
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CRITICAL STRAIN 


During the test the temperature level in the furnace 
was maintained within +2! .°F. by an instrument that 
recorded as well as controlled the temperature level 
utilizing a sensing thermocouple attached to the center 
of the specimen. In all tests, load was applied after the 
temperature level reached equilibrium at 650°F. and re- 
mained so for 1 hour. 


Basic Creep Properties 


Attempts have been made by investigators to cor- 
relate creep buckling theory with experiment using ten- 
sile creep data in place of compressive creep data on the 
assumption that the tensile and compressive creep be- 
havior is identical. However, recent tests? have shown 
that no such reciprocity exists as a general rule, and, 
therefore, it is necessary to determine both types of data 
to prove the applicability of the tensile creep data. 

In accumulating basic creep data, a statistical ap- 
proach is required because of the inherent variability of 
the creep phenomenon. Specimens machined from the 
same material and apparently undergoing that same 
mechanical loading and thermal environment give re- 
sults which vary. As in fatigue studies, therefore, the 
evaluation of such data requires the use of statistical 
methods. Estimates are required of the inherent varia- 
tion in the basic creep properties to permit valid con- 
clusions to be drawn regarding the degree of correlation 
with creep buckling theory. 

The basic creep properties of 3003-0 at 650°F. under 
both compressive and torsional loading were obtained 
under the stress levels used for the creep buckling study 
reported in Section (3). Both 0.060-in. and 0.040-in. 
thick specimens were used to establish the creep proper- 
ties. Although the specimens were permitted to creep 
to buckling. only data from the steady-state region, 
which is assumed to be unaffected by buckling, were 
used. 

An experimental check of potential deformation in 
the threaded portion of the torsion specimens was de- 
termined from a short specimen with 1 2-in. threaded 
portions but with a 1 S-in. gage length of 0.060-in. wall. 
The specimen was loaded to the maximum stress en- 
countered during the program and allowed to creep to 
15° at 650°F. By noting the appearance of the ex- 
trusion markings on the inside of the tube, it was found 
that the scratches in the region of the threads remained 
straight and undistorted while those in the reduced por- 
tion of the specimen were distorted 15°F. The absence 
of any distortion in the threaded portion clearly proved 
the absence of creep in this portion of the specimen. 

The creep properties obtained are summarized in 
Figs. 3and 4. Illustrated is the arithmetic mean of the 
creep data at each stress level. Estimates for statistical 
limits above and below the mean are the +3 standard 
deviations which imply that 99.7 per cent of the data lie 
within these limits for a normal distribution. 


(3) CREEP BUCKLING EXPERIMENTS ON CYLINDERS 


Utilizing the techniques and equipment described in 
the previous section, a series of creep buckling experi- 
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ments were conducted on 3003-0 aluminum tubes at 
650°F. under compressive and torsional loads. The 
pertinent dimensions and stress ranges are summarized 
in Tables | and 2. 

The results of the tests are presented in terms of the 
critical time and critical strain for creep buckling at a 
specified stress level. These data were determined from 
creep curves obtained during the course of the creep 
buckling test and by utilizing the procedure depicted in 
Fig. 5. Typical compressive and torsional creep buck- 
ling specimens before and after testing are illustrated in 


Figs. 6 and 7. 
(4) CORRELATION WITH CRITICAL STRAIN APPROACH 


Compression 


In order to correlate the critical strain approach with 
test data, it is necessary to determine the critical strain 
as discussed in reference 1. It is well known that for 
circular cylinders under axial compressive loads, the re- 
sults of classical, small deflection, elastic stability the- 
ory®* ‘are not in agreement with test data for cylinders 
with Rk h values greater than approximately 500. 
However, for cylinders of relatively small RK h which 
generally are subject to plastic instability, the agree- 
ment between the results of small deflection theory 
modified so as to include plasticity and test data is rela- 


tively good.> Consequently, for the cylinders used in 
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Fic. 3. Compressive creep properties of 3003-0 aluminum at 
650°F. 
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TABLE 1 
Compressive Creep Buckling Specimens 
3003-0 Aluminum Tubes 1.045 in. Inner Diameter at 650°F. 
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h R/h L?/Rh Oa 
0.060 in. 1.0 in. 9.7 28.6 1,650 psi 100.8 hrs. 
0.060 100.1 
0.040 2.0 14.1 177 64.4 
0.040 63.5 
040 63.2 
0.030 2.0 18.4 241 42.1 
0.060 1.0 9.7 28.6 1,925 13.6 
0.060 13.7 
0.040 2.0 14.1 177 8.6 
0.040 8.6 
0.040 8.0 
0.030 2.0 18.4 241 5.9 
0.060 1.0 9.7 28.6 2,185 5.9 
0.060 6.2 
0.040 2.0 14.1 177 + 
0.040 
0.040 2.9 
0.030 £0 18.4 241 2.4 

TABLE 2 
Torsional Creep Buckling Specimens 
3003-0 Aluminum Tubes 1.045 in. Inner Diameter at 650° F. 

h R/h L/R Ta 
0.060 in. 5 in. es 8.6 875 psi 31.4 hrs. 
0.060 
0.040 14.1 8.9 17.3 
0.040 7.2 
0.040 16.8 
0.030 18,4 9.0 12.6 
0.060 5 9.7 8.6 960 1s.1 
0.060 13.3 
0.040 14.1 8.9 7.5 
0.040 
0.040 7.4 
0.030 18.4 9.0 5.5 
0.060 5 9.7 8.6 1,090 6.5 
0.060 6.6 
0.040 14.1 8.9 3.0 
0.040 3.2 
0.040 3.2 
0.030 18.4 9.0 3.1 

this investigation which were in the KX // range of 10 to Torsion 


20 it can be assumed with reasonable confidence that 
the results of small deflection theory may be used. 

From the small deflection theory of elastic buckling, 
the critical strain of long cylinders (L?/ Rh > 9) under 
axial compression is given by 


= (1) 


Since the cylinder is long, this relation applies for all 
boundary conditions between simple support and 
clamped edges. It can be observed from Table | that 
the tubes used in this investigation were in the long 
cylinder category. 

The critical compressive strain has been computed 
for the three groups of cylinders listed in Table 1 ac- 
cording to Eq. (1). The results are given in Table 3. 
In the last column of Table 3 the average experimental 
critical strains are given as obtained from Table |. It 
can be observed that the theoretical and experimental 
values are in good agreement with a maximum dis- 
crepancy of approximately + per cent. 


The elastic buckling 
length under torsional loads is given by the following 


formula :* 4 


€cr 
0.0630 in. 
0. 0605 
0.0419 
0.0416 
0.0420 
0.0313 
0.0608 
0.0600 
0.0410 
0.0415 
0.0412 
0.0325 
0.0605 
0.0615 
0.0425 
0.0395 
0.0405 
0.0306 


Yer 
0.0545 in./in. 
0.0570 
0.0310 
0.0310 
0.0312 
0.0230 
0.0530 
0.0540 
0.03820 
0.0315 
0.0320 
0.0235 
0.0545 
0.0550 
0.0315 
0.0340 
0.0325 
0.0232 


ter = [ew’E/12(1 — RY? URL)? 


for 


where é 


10(h/R)'2< L/R< 3(R/h)'? 


= 0.85 
0.95 


II 


simply supported edges 
clamped edges 


( 


stress of cylinders of moderate 


9) 


The cylinders used in this investigation as listed in 


TABLE 3 


Critical Compressive Strains 


€cr, Eq. (1) 
0.0624 in. /in. 
0.0430 
0.0329 


€cr 
average 
experimental 
0.0611 in. /in. 
0.04138 
0.0315 
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TABLE 4 
Critical Shear Strains 


Yer 


average 
R/h L/R Yer, Eq. (5) experimental 
9.7 8.6 0.0535 in. /in. (). 0547 in. /in. 
14.1 8.9 0.0824 0.0319 
9.0 0.0237 0.0232 


Table 2 fall in the moderate length category. The pre- 
cise boundary conditions in the region between the re- 
duced section and the threaded ends are difficult to 
evaluate. However, in view of the small difference in 
the value of c existing between the limiting conditions of 
simply supported and clamped edges, the following 
average value is used: 

c = 0.89 (3) 


In order to use Eq. (2) for creep buckling, it is nec- 
essary to determine the critical shear strain. For this 
purpose, the following relation is utilized: 

y=7/G (4) 
where G = E/3 


Because of the fact that y represents the shear creep 
strain at buckling (it corresponds to the axial deforma- 
tion in compression), it is assumed that the fully plastic 


value of vy = 0.5 applies, and, therefore, Eq. (2) in 
terms of critical shear strain becomes 

5/4 1/2 

Yer = 2.32 (R/L)' (5) 


The critical shear strain has been computed for the 
three groups of cylinders listed in Table 2 according to 
Eq. (5). The results are given in Table 4. 

In the last column of Table 4, the average experimental 
critical strains are given as obtained from Table 2. The 
theoretical and experimental values are within 2 per 
cent of each other. 


Creep Buckling 

The data presented in Tables 3 and 4 provide a partial 
check of the theory in that the critical strains are pre- 
dicted within acceptable accuracy. However, an over- 
all check of the theory is provided from applied stress - 


STRAIN 
“cr 
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TIME ler 
Fic. 5. Method of determining critical time and critical strain. 
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Fic. 6. Compressive creep specimens. 


Torsional creep specimens. 


Fic. 7. 


critical time presentations of the test data such as 
shown in Figs. Sand 9. 

The test data in these Figures are those given in 
Tables 1 and 2. The theoretical values were estab- 
lished by plotting the theoretical critical strains of 
Tables 3 and 4 upon the basic creep data presented in 
Figs. 3 and 4 and then determining the critical time for 
each stress level. In Figs. S and 9 the heavy line corre- 
sponds to the mean basic creep curve while the outer lines 
represent the 99.7 per cent limits. 

It can be observed from Figs. S and 9 that the pre- 
dictions of the critical strain approach are in relatively 
good agreement with the test data. In most cases the 
test data are within the limits of scatter of the basic 
creep data and generally lie close to the mean curve. 


(5) CREEP BUCKLING UNDER CycLic LOADING 


In piloted aircraft designed for thermal flight, it 
would hardly be reasonable to expect that aerodynamic 
heating will be of sufficient duration for the structural 
element to creep buckle under a single constant load ap- 


plication. A more reasonable situation would be one in 


which the element is exposed to various stress levels and 
elevated temperatures for various periods of time.* 
Under such conditions creep buckling may be an ac- 
ceptable criterion for estimating the life expectancy of 
the aircraft. 
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Fic. 8. Comparison between compressive test data and critical 
strain approach. 


In order to obtain a preliminary insight into creep 
buckling behavior under cyclic loading conditions, the 
following exploratory investigation was conducted. 
A compression and a_ torsion specimen of 3005-0 
aluminum with dimensions given in Table 5 was sub- 
jected to a one-hour-on and one-hour-off loading cycle 
until creep buckling occurred. In determining the crit- 
ical time given in Table 5, only the net time under load 
was totaled. 

A comparison of the critical time (9.4 hours) and 
strain (0.0430) for the compressive specimens under 
cyclic load with the comparable data from the constant 
load tests of Table 1 (8.4 hours) and (0.0412) reveals 
that the cumulative creep hypothesis’ is somewhat 
conservative in this case, although not unreasonably so. 
For torsion, the comparison of cyclic load and constant 
load data are 9.1 hours vs. 7.4 hours and 0.0375 vs. 
0.0318, respectively. Again some conservatism of the 
cumulative creep hypothesis is indicated although the 
data are too few to permit any definitive conclusions to 
be drawn. 


(6) CREEP BUCKLING AND CRIPPLING 
OF FLAT PLATES 


For columns under compression and cylinders under 
compression or torsion, which have been the cases con- 
sidered previously herein, it has been well established 
from short-time room temperature tests that buckling 
and failure are essentially coincident. This is primarily 
due to the fact that no significant tensile membrane 


TABLE 5 
Creep Buckling Under Cyclic Loading 


h R/h 
0.040 in. 2in. compression 
oa = 1,925 psi 


0.040 14.1 5 torsion 9.1 0.0375 
Ta = 960 psi 


Loading net for OF Ver 
9.4 hrs. 0.0480 in./in. 


1100 1100 af 


R/n=18 4 R/h=14 


psi si 
| 
| 
5 10 15 5 10 5 20 
fer Ars ler-hrs 
1000 
psi 
900 
800 
(4) 5 10 1§ 20 25 30 35 


Fic. 9. Comparison between torsional test data and critical 
strain approach. 


stresses arise in the post-buckled state of such structural 
elements. 

For flat plates, however, significant tensile membrane 
stresses can arise after buckling with the result that the 
maximum average strength can appreciably exceed the 
buckling stress if buckling occurs well within the elastic 
range. If buckling occurs above the proportional limit 
then again the buckling stress and maximum average 
strength are essentially coincident. Thus, in con- 
sidering creep of flat plates, it becomes necessary to dis- 
tinguish between creep buckling and creep crippling in 
certain cases. 

In order to obtain some preliminary data on flat 
plates, two square tubes of extruded 3003-0 aluminum 
were subjected to a compressive stress of 2,185 psi at 
650°F. using the equipment described in Section (2). 
The specimens were 4 in. long and 1.25 in. wide, with 
0.078 in. walls. The experimental critical strain for 
creep failure of the two specimens was 0.0153 in. in. 
and 0.0154 in. in. 

Since a square tube is structurally equivalent to four 
simply supported plates, the theoretical critical strain 
was computed from 


= 3.62 (6) 


For the square tube, 6/4 = 15, and therefore the theo- 
retical critical strain is 0.0156 in. in., which compares 
very favorably with the two experimental values. 

The fact that the theoretical critical strain for creep 
buckling and the experimental critical strains for creep 
crippling are essentially identical is attributed to the 
relatively small value of b / of the flat plates used in 
this experiment. This value of } h is sufficiently small 
so that short-time buckling would be inelastic and there- 
fore essentially coincident with the crippling strength. 

The basic compressive creep properties of the square 
tube were not determined and therefore it is not possible 

(Continued on page 458) 
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Effects of a Time-Varying Test Environment 
on the Evaluation of Dynamic Stability With 
Application to Flutter Testing 


WILMER H. REED, III* 
Langley Aeronautical Laboratory, NACA 


SUMMARY 


The paper considers some factors which serve to make the ob- 
served stability of a dynamic system in a time-varying test en- 
vironment differ from its stability under steady testing condi- 
The dynamic response of a system with a time-varying 
and various frequency- 


tions 
excitation frequency is briefly discussed, 
sweep techniques are considered. An analysis of the transient 
motion of a single degree of freedom system whose damping and 
natural frequency change with time indicates that the apparent 
stability depends upon which vibration measurement is being 
observed. Vibration time histories of acceleration are found to 
ye about three times as sensitive to parameter variations as are 
time histories of velocity and displacement. 

The results of the analysis are applied to the problem of flutter 
testing in a time-varying test environment. An approximate 
relation is derived for the magnitude of the apparent shift in the 
flutter boundary due to the unsteady nature of the tests The 
analysis is supplemented by analog-computer studies showing the 
behavior of a single degree of freedom pitching airfoil in a simu- 
lated rocket-sled run and a two-dimensional bending-torsion flutter 
condition in which the torsional stiffness of the airfoil varies with 


time 


SYMBOLS 


b = one-half chord 

¢ = viscous damping coefficient 

Af = change in frequency per cycle 

h = vertical displacement of axis of rotation 


J_ = moment of inertia of airfoil about elastic axis per unit 


span 

J, = Bessel’s function of order n 

K = parameter being varied with time 

k = spring constant 

m = mass 

Jf = Mach Number 

N = ratio of decrement of variable-coefficient solution to de- 
crement of steady-state solution [see Eq. (20) | 

Q = function defined in Eq. (10) 

ro = radius of gyration of airfoil referred to axis of rotation 

= time 

Y= period of oscillatory motion 

Vo = velocity of airfoil 

x = displacement of mass in Eq. (5) 

x9 = dimensionless distance in chord lengths from leading 


edge to axis of rotation 
Xq@ = dimensionless static unbalance in semichord lengths 
about axis of rotation (positive for rearward center of 


gravity) 


¥ = variable defined by Eq. (8) 
@ = angular displacement of axis of rotation 
@ = nondimensional frequency-sweep parameter 


[see Eq. (1) ] 


Presented at the <Aeroelasticity-II Session, Twenty-Sixth 
Annual Meeting, IAS, New York, January 27-30, 1958. 
* Aeronautical Research Engineer. 


amplitude at ¢ = 0 


= decrement, 
amplitude att = 7 

¢ = damping ratio in per cent of critical damping 

AX = rate-of-change-of-frequency parameter 

a = relative density of typical section, m/mpb* 

p = air density 

w = angular frequency 

w, = angular frequency of wing in deflection 

We = angular frequency of wing in torsion 

w, = natural angular frequency of single degree of freedom 


vibrating system 


Subscripts 
0 = condition at the beginning of a cycle of oscillation 
= flutter 
vy, é, é = vibration time history of displacement, velocity, and 


acceleration, respectively 


INTRODUCTION 


ECAUSE OF THE COMPLEXITY of modern aircraft 
B structures and uncertainty of the nature of un- 
steady air forces acting on the airplane, it has become 
analytical predictions of 
The degree of 


necessary to supplement 
flutter with some form of flutter tests. 
stability observed in the flutter test is then generally 
expressed in terms of such familiar stability parameters 
as natural frequency and damping ratio in per cent of 
critical damping, the amplitude and phase angle of the 
forced response over a range of frequencies, the power 
input to the exciter, and so forth. It has been found 
that a large variety of flutter problems fall within the 
framework of linear dynamic systems defined by dif- 
ferential equations having constant coefficients. It is 
natural, therefore, to appraise the dynamic charac- 
teristics observed in the flutter tests on the basis of 
methods developed for linear systems with constant 
parameters. The testing technique, however, may be 
one in which some of the system’s parameters vary con- 
tinuously with time during the course of the test. 
This has the advantage of greatly shortening the time 
required to cover a given range of variables but, on 
the other hand, the effect of the rate of change of the 
parameters may completely alter the apparent stability 
of the system as compared with its “‘steady-state”’ 
stability. 

When employing the frequency-response technique 
it is often desirable to let the excitation frequency vary 
slowly with time over the range of frequencies of in- 
terest. The time-varying parameter of importance 
here, assuming that other factors are constant, would be 


the rate of frequency sweep. 


| : 
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435 


436 JOURNAL OF THE AERO/SPACE SCIENCES—JULY, 1958 


The evaluation of stability from transient-response 
data can be altered when unsteady testing conditions 
cause the parameters which define the system to vary 
with time. These parameters may be associated with 
either the test medium or the physical properties of the 
structure being tested. Examples of flutter-testing 
techniques in which the test medium varies with time 
are rocket-sled tests wherein the forward velocity varies 
and blowdown-wind-tunnel tests in which the density 
varies. Structural variations with time may be due to 
such factors as rapid transient heating which would 
cause the stiffness to change, or time variations in mass 
or center of gravity position. 

At the present time, the flutter engineer must rely 
on previous experience or trial and error to ensure that 
the rates of change of parameters which affect flutter 
are not large enough to cause appreciable differences 
between the observed flutter boundary and the “‘steady- 
state’ boundary. Laidlaw and Beals have commented 
in reference | that early rocket-sled tests showed that 
rapid accelerations into previously unexplored speed 
ranges are likely to yield misleading results. 

It is the purpose of this paper to study how and to 
what extent apparent stability can be affected by un- 
steady testing conditions. On the basis of some pre- 
vious work relating to effect of frequency sweep on the 
response of a dynamic system, a sweep technique is sug- 
gested which makes the error due to sweep independent 
of frequency at which resonance occurs. The effects 
of time-varying parameters on the apparent stability 
measured from the transient motion are analyzed for a 
single degree of freedom system, and the results of the 
analysis are used to predict shifts in flutter boundaries 
and changes in modal damping due to unsteady testing 
conditions. Illustrations are given for the transient 
motion of a single degree of freedom pitching airfoil 
undergoing forward acceleration, and the bending-tor- 
sion flutter of an airfoil with time-varying torsional 
stiffness. 


ANALYSIS 

Effects of Variable Excitation Frequency 

The frequency-sweep technique, in which the excita- 
tion frequency changes continuously with time, is 
widely used in airplane flight flutter and dynamic-sta- 
bility tests. When employing this technique it has 
been observed that the maximum amplitude of re- 
sponse is less than the response that would occur if the 
frequency were held constant on resonance; also, the 
apparent resonant frequency is shifted, being less than 
the true resonant frequency when the frequency is de- 
creasing with time and greater when the frequency is 
increasing. These effects are illustrated in Fig. 1 for a 
single degree vibrating system in which the excitation 
frequency increases and decreases linearly with time. 

The subject has been treated by several authors in 
the fields of electrical and mechanical engineering. The 
first known analysis relating to this subject was by 
Lewis,? who treated the problem of accelerating rotat- 


FREQUENCY 


FREQUENCY 


RESPONSE 
4- _—-STEADY STATE 
2r 
1 1 1 1 j 
8 9 10 Ll 12 13 
FREQUENCY 


Fic. 1. Effect of variable excitation frequency on the response 


of a single degree of freedom vibrating system. 


ing machinery through a resonant frequency. Hok in 
reference 3 made a similar analysis with application to 
electrical networks and extended Lewis’ results by in- 
cluding the effects of frequency sweep on the phase 
angle. 

A particularly convenient form of the solution for a 
lightly damped single degree of freedom resonance sys- 
tem excited by a linearly varying frequency is pre- 
sented in reference 3. Here it is found that differences 
between the steady-state response and the response 
when the excitation frequency varies with time is re- 
lated to a single nondimensional parameter defined as 


& = fw,/V/ 2a (1) 


where ¢ is damping ratio in per cent of critical damping, 
w, 1S the undamped natural frequency, and w is the rate 
of change of exciter frequency. The larger the fre- 
quency-sweep parameter &, the smaller will be the 
errors due to frequency sweep. Data in reference 3 
indicate that reasonable values of &@ for dynamic-sta- 
bility tests fall between 1.4 and 4.0, giving errors in the 
apparent peak amplitude response between about 10 
per cent and less than 1 per cent, respectively. Similar 
errors are indicated for the apparent resonant frequency. 

If several modes are present in a given range of fre- 
quencies, the & criterion for frequency sweep suggests 
the possibility of varying the frequency in a manner 
such that & will always be equal to or greater than an 
assigned minimum value regardless of the resonant 
frequencies of the modes being studied. Thus, with an 
estimate of the minimum damping likely to be en- 
countered in the tests, the variation required to main- 
tain a constant & may be derived from the definition 
of & in Eq. (1): 


dw(t)/dt = (1/2)d w(t) (2a) 


where \ = (€nin/®)? and (min iS the minimum damping. 
(Eq. 2) shows that \ is proportional to the per cent of 
change in frequency per cycle: 

A = (1/m) (Af/f) (2b) 


Integrating Eq. (2a) and letting the exciter frequency 
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EFFECTS OF A 


at / = 0 be we gives a variation of frequency with time 
such that the per cent of change in frequency per cycle 
is constant.* 

w(t) = aw {1 — (1 2) wo At] (3) 


It is interesting to compare Eq. (3) with other pro- 
gramed frequency variations which are often used in the 
frequency-sweep technique. Fig. 2 shows an example 
of how a frequency range of, say, from 5 to 35 eps 
might be surveyed using a linear, a logarithmic, and a 
constant per cent frequency change per cycle sweep 
technique. It is assumed that at ¢ = 0 the frequency 
and rate of change of frequency are the same for the 
three methods shown. The variation of frequency with 
time is given by 


w(t) = + (1/2) (4a) 


for the linear sweep, by 


(1/2) Agwot 
w(t) = ave 


(4b) 


for the logarithmic sweep, and by Eq. (3) for the con- 
stant per cent of frequency change per cycle (A = 
constant) sweep relation. 

The conditions assumed are & = 2.0 and ¢,,;, 
0.05, giving \ = 0.000625 or Af/f ~ 0.002. The most 
obvious point to make about the Figure is the time re- 
quired to sweep from 5 to 35 eps by the different tech- 
niques shown. The times are 1.5 min. for the curve of 
\ = constant [Eq. (3)|, 3.2 min. for the logarithmic 
sweep, and 10 min. for the linear sweep. It should be 
noted that as the frequency range to be surveyed is re- 
duced, the sweep times required by the three methods 
shown become more nearly equal. 

If the sweep parameter & is evaluated on the basis 
of the instantaneous values of w(f) and @(/) during the 
frequency sweep, it is found that & increases as V w(t) 
for the logarithmic sweep and increases as w(t) for the 
linear sweep. This would indicate that the higher the 
exciter frequency the smaller the errors due to sweep 
for the linear and logarithmic frequency variations; 
however, if & is selected to give acceptable errors at the 
low frequencies there is no need to improve the accuracy 
at higher frequencies, especially since the cost of better 
accuracy is the increased time required to complete the 


tests. 


Transient Response With Time-Varying Parameters 


Recent studies of the dynamic behavior of missiles 
and airplanes undergoing rapid altitude or speed 
changes have shown the importance of taking into ac- 
count the time rate of change of parameters affecting 
the system. For example, Zimmerman in reference 4 
shows that in ascending flight a disturbance may result 
in oscillations which grow in amplitude sufficiently to 
cause the vehicle to leave the atmosphere in a tumbling 


“It should be pointed out that for systems having more than 
one degree of freedom and for frequency variations other than 
linear, the errors indicated by the @ criterion in reference 3 are 


only approximate. 


TIME-VARYI 


NG ENVIRONMENT 437 


TEST 


motion even though its steady-state dynamic charac- 
teristics are known to be stable. Conversely, Friedrich 
and Dore indicate in reference 5 that in re-entering the 
atmosphere the maximum angle of attack of a missile 
is reduced even when aerodynamic damping is ne- 
glected. 

These predictions stem from the fact that the equa- 
tions of motion defining the system contain coefficients 
that vary with time. In missile applications the effects 
of time-varying coefficients are of importance in predict- 
ing the maximum loads, temperatures, trajectories, etc. 
In testing for dynamic stability, however, the effects 
of a time-varying test environment are important be- 
cause the interpretation of stability may be affected 
by the unsteady nature of the tests. 

To evaluate the effect of time-varying parameters on 
the apparent stability of a dynamic system, a study is 
made of the transient motion of a single degree of free- 
dom system with time-varying damping and frequency. 
The analysis somewhat parallels that of reference 5; 
however, in this paper the primary emphasis has been 
placed on the apparent damping of the motion. 

Equation of Motion—Consider the familiar damped 
single degree of freedom vibrating system defined by 
the equation of motion, 


+ + = O (5) 


where m/(t) is the mass, c(t) a viscous damping coeffi- 
cient, k(t) a spring constant, and x the deflection of the 
mass from its equilibrium position. The coefficients 
are assumed to be gradually varying functions of time 
such that over any given cycle of motion the variation 
is approximately linear. The coefficients in Eq. (5) 
then become 
m(t) = mo + mola 


c(t) = Col (0) 
R(t) ko + Rot ) 


I 


where the 0 subscripts denote conditions at the begin- 
ning of the cycle. The time interval of interest in the 
analysis to follow will be of the order of the period of 
oscillation. Incremental values of the coefficients in 
Eq. (6) in one vibration cycle are assumed to be so small 
that products between the increments are negligible. 
Dividing Eq. (5) by the m(t) defined in Eq. (6) and drop- 
ping the second-order terms gives 
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KA + + wo? [1 + 2 lx = O (7) The second transformation introduces Q [Eq. (11)| 
as a new independent variable which, when substituted 
where : 
into Eq. (9), gives 
Co = wy? = ko/mo, Me 
(d*t/dQ*) + [Q/(2a@)*?] # = O (12 
(o/'Ro Eq. (12) may be identified as a form of Bessel's dif. 
Solution—-Two transformations of variables are use- ferential equation (see reference 7), having as its solu. 


ful in obtaining a solution to Eq. (7). The first trans- 
forms the dependent variable’ by means of the relation 


= ,—w0[to-+ (fot/2 


which, when substituted into Eq. (7), eliminates the ¥ 
term to give 

x + Q(t)x = 0 (9) 
where 


Q(t) = wo?) {1 wo) | 


2 [ (wo a) |t (10) 


Assuming the system to be lightly damped (¢) << 1) 
and again neglecting the second-order €) and a) terms 


tion 


BJ (13 


where J/..;/3 is Bessel’s function of the first kind of 
order 1/3 and 4 and Bare arbitrary constants. Eq, 
(13) can be handled more conveniently by replacing 
/,(y) by its asymptotic expansion for large values of the 
argument y. 
like a damped trigonometric function: 


The approximation for J,(y) behaves 


cos[y — (2/4) — (mm 2)] 
Jily) = (14 
2)y 
The necessary requirement that the argument 


Q* *3wu, in Eq. (13) be large is consistent with the 


gives 
; assumption that the variation in coefficients per cycle 
Y(t) wo? [1 — (Fo/wo) 2a/wy)t (11) 
wo" | is small. 
When Eq. (14) is used, the working form of Eq. (13 
in terms of the original variables x and ¢ becomes 
xv A + (ot cos } (wo? 3a) [1 = (Fo ) ¢} 
x(t) = - (15) 
(1 (fj ao) + (wo wy) t]' 
where .! and ¢ are arbitrary constants. 
Period and Envelope of Motion—The motion of the : 
mass as given by Eq. (15) is an oscillation of vary- follows: 
ing frequency and amplitude. An instantaneous fre- log A = 2x¢ (19 
quency may be defined by differentiating the argument . 
of the cosine term to give amplitude at ¢ = 0 
where A= = 


w(t) = — (Fo/ao) + (2ce/w)t (16) 
The period of motion 7 may be obtained by finding 
the time required to increase the argument by 27, 
which, to a first-order approximation, is 


= (27 wo) {1 = (a wo”) + (fo 2wo) | (17) 


As the phase angle ¢ in Eq. (15) is varied, a family 
of oscillating curves is generated. The envelope con- 
taining the curves is 
[x(t)| = |x(O)| [1 + (18a) 


where 1x(0)| is a point on the envelope at time ¢ = 0. 
Similarly, envelopes for |x(t)| and |;¢(t) are found from 
Eq. (15) to be 


= [1 + + gp) 


Il 


-1pparent Damping—Knowing the period and enve- 
lope of the motion, it is possible to compute the loga- 
rithmic decrement in a manner analogous to that for con- 
stant-coefficient systems. For lightly damped systems, 
the logarithmic decrement is related to the damping as 


amplitude at ¢ = 7 


An important difference should be noted here be- 
tween the damping of the constant-coefficient system 
and the apparent damping of the variable-coefficient 
system. For the constant-coefficient system, time his- 
tories of x and all derivatives of x reveal the same 
damping. On the other hand, when discussing the 
damping of the variable-coefficient system it is neces- 
sary to specify the particular variable being considered. 
This can be seen from Eqs. (18) by noting that the 
ratios 


x(t)!’ 


x(0)| 
ix(t)|’ 


| 


are different in each case. It is possible, then, for the 
time history of x(¢), say, to appear damped while over 
the same time interval X(f) or ¥(f) may be diverging. 
This point is illustrated in Fig. 3, which shows the time 
histories of x, x, and * for an undamped single degree 
of freedom system whose natural frequency is increas- 
ing at the rate of a 10 per cent change in frequency per 
cycle. The effect of increasing frequency is to make the 
undamped system behave as both a stable and an un- 


stab 
tory 
peal 
hist 
bot! 
fine 
bee! 
mat 
1 
solt 
div! 
vel 
“st 
me! 
tha 
are 
the 
ste: 
ter 
ren 
for 
Th 


E 
# 
the 
Th 
dis 
vel 
ac 
wl 
th 
| 
It 
° 
in 
et 
3 
| 
it 
d 


4- (11)] 
stituted 


(12 


dif- 
Soly- 


(13 
ind of 
Eq. 
lacing 
of the 
‘haves 


(14 


iment 
n the 
cycle 


(19) 


EFFECTS OF A TIME-VARYING TEST ENVIRONMENT $30 


stable system, depending upon the particular time his- 
tory being observed. The displacement oscillation ap- 
ears stable, whereas the velocity and acceleration time 


This simultaneous existence of 


histories are diverging. 
both stable and unstable characteristics of a system de- 
fined by equations with time-varying coefficients has 
been noted by Grensted in reference 8, in which various 
mathematical definitions of stability are discussed. 

The apparent damping of the variable-coefficient 
solutions may be expressed in a convenient form by 
dividing the decrements evaluated from displacement, 
velocity, and acceleration solutions [Eqs. 18)| by a 
“steady-state” decrement. The steady-state decre- 
ment is defined as the decrement obtained by assuming 
that the coefficients in the equation of motion [Eq. (7) | 
are constant and equal to their average values during 
the cycle in question. With this definition for the 
steady-state decrement, it is found that the exponential 
terms are identical for both the variable-coefficient dec- 
rements and the steady-state decrements and, there- 
fore, cancel when forming the ratio between the two. 
The resulting decrement ratios 

Vv = A variable coeflicient 
A steady state) 

then become simple functions of a single parameter X. 
These relations for three vibration measurements are 


displacement : NV, = (1 + 2nd)" (20a) 
velocity : N,; = (1 + (20b) 
acceleration : Ny = (1 + (20c) 


where the subscripts denote the time history from which 
the decrement is derived and 


A = 2a wo” = (1 1) (Af f) (21) 


It is interesting to note that A, or the per cent of change 
in frequency per cycle, was also the significant param- 
eter in determining the frequency-sweep rate for 
forced oscillations [see Eq. (2) ]. 

From Eqs. (19) and (20) the apparent damping is 


= (1/27) log N, state) 
OISPLACEMENT VELOCITY ACCELERATION 
| 
a 
| 
y + | 
| 
TIME 


_ Fic. 3. Time histories of an undamped system—effects of 
increasing natural frequency with time on the motion of an un- 
damped single degree of freedom system. 
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Fic. 4. Damping increment due to time-varying frequency 
or (1/2) log N, = fn — [fo + (For, wo) (22 


where 7 is x, X, or ¥, depending on the time history used, 
and ¢) + (for wy) is the damping ratio averaged over 
the cycle. The term (1/27) log .\, is, therefore, the 
difference between the apparent damping and _ the 
averaged steady-state damping. The damping-ratio 
increment (1 27) log .V, is plotted against the per cent 
of change in frequency per cycle in Fig. 4. To illus- 
trate the Figure, consider the time histories for the un- 
damped system shown previously in Fig. 3. Here the 
frequency is increasing at the rate of 10 per cent per 
cycle; therefore, from Fig. 4 the apparent damping of 
displacement, velocity, and acceleration is ¢, = 0.007, 
¢; = —0.007, and ¢; = —0.02. The damping incre- 
ments for velocity and displacement are approximately 
equal and opposite and are about one-third as large as 
the increment for acceleration. As compared with the 
steady-state damping, vibration measurements of veloc- 
ity and acceleration would appear less damped for in- 
creasing frequency and more damped for decreasing fre- 


quency. The opposite would be true regarding the dis- 
placement time histories. 
Illustrative Example—The following example will 


serve to illustrate and summarize the foregoing analysis. 
Consider a single degree of freedom dynamic system de- 
fined by a constant-coefficient differential equation. 
Assume that the characteristic damping and frequency 
of the system are known functions of some parameter 
K. If K is varied in a prescribed manner A(t), and the 
free vibrations of the system are observed, in general 
there will be a difference between the apparent damp- 
ing and the damping which would exist if conditions 
were steady at each instantaneous value of A(t). To 
estimate how much the observed damping would be af- 
fected by varying A with time, the following procedure 
is outlined. 

(1) Plot ¢ and w against A. 

(2) Graphically or analytically determine the slopes 
(O¢ OA)y and (Ow/O0A), for various values of Ay. 

(3) From the slopes and the prescribed rate of change 
of A with time, compute 


Gy = (Ow OK), Ay, 


fo (O¢ OK Ko, A= 2a) 
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Fic.. 5. 


Modal characteristics of a single degree of freedom pitching airfoil as functions of reduced speed parameter. 


mrpbt/Ia = 0.0016; pitch axis at leading edge. 


(4) Compute the average steady-state damping ratio 
and the average A; 


Sar = (for wy), Ko + (Ar we) 


(5) Using results obtained in steps (3) and (4) above, 
compute apparent damping for displacement, velocity, 
and acceleration from Eq. (22) and plot against A,,. 


APPLICATIONS OF ANALYSIS TO FLUTTER TESTS 


The disturbed motion of an elastic structure in an air 
stream is governed by the simultaneous action of elas- 
tic, inertia, and aerodynamic forces. These forces are 
a result of the motion itself and, for small disturbances, 
are linearly dependent on the motion. Linear sta- 
bility theory shows that the free motion of a dynamic 
system may be regarded as the superposition of natural 
modes, each mode having associated with it a par- 
ticular damping, natural frequency, and deformation 
shape. Thus, the disturbed motion of a dynamic sys- 
tem may be completely described in terms of the 
natural modes of the system. 

When considering flutter, the natural modes stem 
directly from the structural ground-vibration modes. 
If the air speed is increased from zero these modes are 
altered by the action of aerodynamic forces. Even- 
tually, a speed may be reached where the damping of 
one of the modes vanishes. The disturbed motion of 
the system would then be a self-sustained oscillation at 
the frequency of the undamped mode. The air speed 
for which this neutrally stable condition exists is called 
the flutter speed, and the natural mode involved is the 
flutter mode. 

In order to estimate the effects of time-varying test 
conditions with regard to flutter testing, it will be as- 
sumed that the normal vibration modes may be treated 
individually as single degree of freedom systems. In 

other words, a normal mode, whose damping and fre- 


quency are related to some physical parameter in a 
known manner, will be replaced by a single degree oj 
freedom system having identical frequency and damp- 
ing characteristics. 


Apparent Shift in Flutter Boundaries 


The previous analysis has shown that if a condition 
of neutral dynamic stability is approached while some 
parameter affecting the stability is continuously varied, 
the apparent stability boundary may differ from the 
actual boundary that would exist if conditions were 
steady. Further, the apparent boundary would de- 
pend upon whether stability was judged from a vibra- 
tion measurement of displacement, velocity, or accelera- 
tion. 

To obtain an approximate formula for the magnitude 
of the shift in flutter boundaries due to a time-varying 
test condition, consider the following example. As- 
sume that the variation of damping and frequency with 
some parameter A in the vicinity of flutter can be ap- 
proximated by a linear relation, such that 


OK ), = (O¢ OK), (23 


wo = Wy + (Ow 


where AA, is the increment in A from its value at the 
flutter boundary and the subscript f refers to condi- 
tions at the flutter point. Further, assume that A 
varies linearly with time through its critical value. The 
rate of change of frequency and damping with time is 
then 

ty = (0w/dK), AK, (24 


It can be seen from Eq. (22) that the condition for ap- 
parent neutral stability (¢, = 0) is satisfied when 
(1/27) log V, = — [fo + (For/wo) (25) 


Substituting Eqs. (20), 
gives an expression for 


(23), and (24) into Eq. (25) 
the shift in flutter boundary 
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AK., = Ake + | rh w,) for each of the three vibration 
measurements of interest. An explicit solution for 
\K,, cannot be obtained; however, by expanding the 
logarithmic term and neglecting higher order terms, the 


following approximate relation is found: 


AAy, = |, (26) 
where, for displacement, p = —1/4; 
» = 1 4; and for acceleration, p = 3 4. 
This brings out the point made previously that ac- 


for velocity, 


celeration measurements are three times as sensitive to 
the effects of time-varying conditions as are the dis- 
placement and velocity measurements. Also, dis- 
placement and velocity time histories would indicate ap- 
parent boundaries that are shifted approximately equal 
amounts on either side of the steady-state boundary. 
This suggests the possibility of predicting the steady- 
state boundary by averaging the boundaries deter- 
mined from velocity and displacement time histories in 
an unsteady test. Finally, it is seen that the larger 
the value of \ and the smaller the slope of the damping 
at the flutter point, the greater will be the spread be- 
tween the actual and the observed flutter boundary. 


Single Degree of Freedom Pitching Airfoil 

Several specific examples of the effects of time-vary- 
ing parameters on the results of flutter tests will now 
be considered. The first illustration is for a single 
degree of freedom pitching airfoil with the axis of rotation 
at the leading edge. While there is some question re- 
garding the practical significance of this mode of flutter, 
the example chosen provides a convenient example for 
showing how the apparent stability in flutter testing 
may be altered by time-varying parameters. Since, in 
general, the disturbed motion of the airfoil is a damped 
oscillation, the aerodynamics of the problem are ex- 
pressed in terms of Wagner's indicial functions for in- 
compressible flow,* making use of the approximations 
for these functions given in reference 10. Variation of 
the steady-state damping and frequency ratios with 
reduced speed coefficient is shown in Fig. 5(a) for a 
density ratio rpb* Ja = 0.0016. As the velocity param- 
eter increases from zero, the damping ratio increases 
to a maximum and then gradually diminishes, passing 
through zero or flutter at |) bwa = 36. The frequency 
ratio remains constant at unity until the point of maxi- 
mum damping, where it begins to increase with |” bwa, 
and it reaches 1.6 at the flutter point.* 

Fig. 5(b) shows the apparent damping that would be 
observed in time histories of displacement, velocity, and 
acceleration when the reduced speed parameter in- 
creases with time at a constant rate d( I" bwa)/d(wat) = 
1.0, where wa is the undamped natural frequency of the 
airfoil about its pivot axis. For 6 = 0.5 ft. and wa = 


“The flutter speed and frequency ratio shown here differ 
from the values shown by Runyan in reference 11. he dif- 
ferences are believed to be due to the approximations used in 
defining the indicial lift function, but for the purpose of illustrat- 
ing the effects of time-varying coefficients, these differences are 


probably unimportant. 
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20 rad. per sec., the forward acceleration would be 
about 6g for the case shown. An acceleration effect 
not considered here is the change in torsional stiffness 
that would occur if the airfoil is not balanced about its 
pivot axis 

Compared with the steady-damping relation, the ap- 
parent damping for increasing speed shows that the 
velocity and acceleration are less damped and the dis- 


placement more heavily damped. The curves for 
variable-speed damping tend to merge together at 
1’ bwa = 10 where the rate of change of frequency be- 
comes zero [see Fig. 5(a)|. This is because differences 
between the three unsteady solutions are due only to 
the rate of change of frequency [see Eq. (22)|. The ap- 
parent damping observed when only the damping coeffi- 
cient varies with time is the average value of damping 
during the cycle and is the same for all time histories 


shown. 


Simulated Rocket-Sled Flutter Tests 


In order to illustrate how the damping characteristics 
predicted in Fig. 5 would appear in the form of time 
histories for a specific flutter test, a rocket-sled flutter 
run for this configuration was simulated on an analog 
computer. The problem was programed on the com- 
puter by the method of reference 12 with minor modifi- 
cations to allow for variations of the velocity with time. 
The results are shown in Fig. 6 as vibration time his- 
tories of the airfoil pitch angle, pitching velocity, and 
pitching acceleration. Also shown is the variation of 
the sled velocity with time. 

In the example shown, a is deflected from equilibrium 
and released at the start of the run. The airfoil is 
rapidly accelerated d( 1° bwa) d(wat) = 1.6 from rest to 
a speed of | bwa = 60, passing through the flutter 
boundary at I) bwa = 36. The speed is then held 
constant for a brief interval, followed by a coast phase 
d(V/bwa)/d(wat) = —0.32 in which the flutter bound- 


SIMULATED ROCKET-SLED FLUTTER TEST 
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Fic. 6. Analog-computer time histories of a single degree of 
freedom pitching airfoil in a simulated rocket-sled flutter test. 
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(b) Damping of torsional mode; torsional stiffness decreasing 
with time = —O0.06. 
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(c) Damping of torsional mode; torsional stiffness increasing 
with time d[wa(t)/wn}?/d(ant) = 0.06. 


Fic. 7. Modal characteristics of bending-torsion airfoil as a 
function of stiffness ratio. WJ = 5.0, uw = 200, V/beo, = 8.0, 
xo = 0.5, Xa = 0.05, ra? =1/4. 
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ary is again crossed in going from an unstable to 4 
stable condition. 
it is interesting to note that through the entire range 


During the high-acceleration phase 


of V/bwa covered, the velocity and acceleration time 
histories diverge, whereas the displacement time history 
is damped. The existence of a stability boundary ot 
\’/bwa = 36 is completely masked by the effects of for. 
ward acceleration. An indication of how the frequency 
is changing over the speed range is shown by the dis. 
Over the 
speed range shown, the frequency is increased by a 


tance between successive oscillation peaks. 


factor of 2.5 

During the second phase, the speed is held constant 
at V/bwa = 60. 
speed and the test conditions are now steady, all time 


Since this velocity is above the flutter 


histories are seen to diverge at the same rate. 
Finally, during the coast phase of the test, the ap- 
parent flutter boundaries indicated by a zero slope of 
the envelope occur at different speeds for the three 
curves shown. As predicted by the previous analysis, 
the apparent boundary judged from the acceleration 
time history is seen to be farthest from the actual 
boundary. 
and velocity measurements are displaced by approxi- 


The apparent boundaries for displacement 


mately equal amounts on either side of the steady 
boundary. 


Bending-Torsion Flutter With Time-Varying Torsional 

Stiffness 

The final example to be considered is the bending- 
torsion flutter of a supersonic airfoil in which the tor- 
sional stiffness varies with time as a result of, say, rapid 
transient heating during the flutter tests. The aero- 
dynamic forces used in the problem are given by piston 
theory'® for a flat plate at a Mach Number of 5.0. 
Other parameters used are listed in the caption of 
Fig. 7(a), in which the steady-state variation of the 
modal characteristics with the stiffness ratio (wa/w;)* 
is shown. Here it is seen that, as the stiffness ratio is 
reduced from an initial value of 4.0, an instability de- 
velops at (wa/w,)” = 1.2 in the mode originating from 
the uncoupled torsion mode w,. 

The apparent damping evaluated from transient 
time histories of the torsion mode is shown in Fig. 
7(b) for (w, w,)* decreasing linearly with time and 
in Fig. 7(c) for (w,/w,)? increasing linearly with 
time. These conditions might approximate a flutter 
test involving elevated temperatures where the tor- 
sional stiffness is first reduced as a result of buckling 
associated with nonuniform heating and then re- 
stored as the structural temperatures become uni- 
form (see reference 14). The rate of 
stiffness with time for the example shown is +0.06 
per rad. based on the uncoupled bending frequency «». 

The points on the Figures are obtained from the 
logarithmic decrement of successive cycles of the analog- 
computer solution for the two degrees of freedom sys- 
tem and are to be compared with the solid lines which 
are predictions of the damping derived from steady- 
state characteristics of the torsional mode. It is seen 
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in Fig. 7(b) that when the torsional stiffness decreases 
with time, the damping measured from the acceleration 
transient is several times greater than the steady-state 
damping. However, when the torsional stiffness in- 
creases with time |Fig. 7(c) |, the acceleration transient 
diverges over the entire range of variations considered. 
\s in the previous examples, Fig. 7 verifies that, for a 
given value of the time varying parameter, in this case 
wa/o the steady-state damping is very closely 
approximated by the mean value of damping obtained 
from transient time histories of angular displacement 


and velocity. 
Other Analog-Computer Flutter Cases 


With regard to the examples chosen to illustrate the 
effects of time-varying parameters on flutter, it should 
be stated that the particular cases shown were selected 
irom a group of analog-computer flutter studies be- 
cause the points to be illustrated were most evident in 
these cases. Many of the analog solutions showed 
only small differences between the steady solution and 
the solution with time-varying parameters because the 
rate-of-change-of-frequency parameter \ was too small 
to have an appreciable effect on the damping. In cases 
where the flutter mode changed abruptly from a well- 
damped to a violent divergent oscillation for a small 
increment variation in the parameter [large (Of OA),, 
see Eq. (26) |, the shift in flutter boundary due to time 
rate of change of the parameter was generally small. 
In other words, on the basis of these results it appears 
unlikely that a violent flutter mode would go un- 
detected or be appreciably altered as a result of un- 


steady testing conditions. 


CONCLUDING REMARKS 


The present paper considers some effects of time- 
varying test conditions on the apparent stability of the 
dynamic system being tested. The analysis is con- 
cerned primarily with the effects of time-varying pa- 
rameters on the transient motions of a system; however, 
the subject of forced response with a time-varying ex- 
citation frequency is briefly discussed. 

The damping assessed from the transient motion of a 
system whose parameters are varying with time is found 
to be most affected by the time rate of change of the 
frequency of the system. Vibration time histories of 
acceleration are shown to be about three times as sensi- 
tive to parameter variations as are the time histories of 
velocity and displacement. Equal and opposite dif- 
ferences are noted between the steady-state damping 
and the damping assessed from velocity and displace- 
ment transient osciliations. 

The results of the analysis are then applied to the 
problem of flutter testing in a time-varying test environ- 
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ment. An approximate relation is derived for the 
magnitude of the apparent shift in flutter boundaries 
due to time variations of any parameters affecting the 
frequency and damping of the flutter mode. It is 
found that the larger the per cent of change in frequency 
per cycle and the smaller the slope of the damping 
curve as the damping vanishes, the greater will be the 
shift in flutter boundary due to time-varying test condi- 
tions. For a given unsteady test environment the 
per cent of change in frequency per cycle will be larger 
for low-frequency than for high-frequency systems; 
consequently, the effect of the time-varying environ- 
ment on stability will be more apparent for the low- 


frequency system. 
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SUMMARY 


The laminar and turbulent mixing of a nonuniform free stream 
with fluid at rest is investigated. Laminar and turbulent vor- 
ticity numbers are defined, which measure the relative impor- 
tance of the free-stream vorticity as compared to the average 
vorticity present in the wake and which increase with the square 
root of x and with x, respectively. Solutions of the subject prob- 
lem are given to a first approximation in the vorticity numbers. 
In comparing the characteristics of the mixing due to a nonuni- 
form stream “; = uy + wy with those due to the uniform stream 
u = up, it is found that vorticity effects over a given length are 
nearly always negligible in the laminar case but became sizable 
in the turbulent case. This is due to the usually higher rate of 
increase with x of the turbulent vorticity number. Results of 
the analysis, as applied to streams having constant vorticity over 
finite y length, indicate the practicability of the use of outer vor- 
ticity as a means to control the rate of transfer of properties be- 


tween two different streams 


SYMBOLS 


= width of the wake 

= nondimensional stream functions 
= upward mass flux per unit length 
= velocity components 

= free-stream velocity = uw + wy 
= Cartesian coordinates 

= nondimensional stream functions [Eq. (18)] 

= eddy coeflicient of kinematic viscosity = eqtox 
= Blasius variable [Eq. (5)] 

= kinematic viscosity 

= density 

= similarity variable [Eq. (15)| 

= turbulent ‘‘vorticity number” |Eq. (17)| 
= stream function 


= laminar ‘‘vorticity number” [Eq. (6)| 


= free-stream vorticity 


INTRODUCTION 


_ sTuby of mixing phenomena has received wide 
attention in recent years. Both theoretical and 
experimental investigations are available to date which 
cover a rather large field of interest.!. There has not 
been, however, any concern in investigating the effects 
of nonuniformity in the free streams upon the char- 
acteristics of mixing. 
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Incompressible Mixing of a Shear Flow 
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Yet, on one hand, practical instances of mixing are 
more likely to take place among nonuniform streams 
than not. 
ing phenomena occurring (a) in internal flow machines. 
(b) as a consequence of the interaction of strong curved 
shock, and (c) in some problems concerning the injec- 
tion of secondary streams into primary, rotational 
streams for cooling purposes. On the other hand. 
there is a variety of cases where the possibility of con- 
trolling the transfer of properties between the two 
streams plays an essential role in the establishment oj 
efficient processes. Such cases arise, for instance, in 
internal flow aerodynamics, in the problems of flame 
stabilization,* flame diffusion* * and combustible miy- 
ture ignition. 

Thus, at least two points appear to be worth in- 
vestigating : 
outer streams changes the characteristics of the miy- 
ing of the corresponding irrotational streams, and (b 
what is the possibility of using vorticity as a means of 
controlling the rate of transfer of properties between 
different streams. 

As a first step toward a fuller investigation of the 
effects of nonuniformity in the free stream, the pres- 


(a) how the presence of vorticity in the 


ent paper deals with the incompressible mixing of a 
nonuniform, constant vorticity, free stream with fluid 
at rest. Both laminar and turbulent cases are con- 
sidered. 

Since the time Ferri and Libby’ first pointed out the 
possible importance of outer stream vorticity on the 
characteristics of dissipative regions several authors 
have been concerned with some aspects of this problem. 
Li® * * has treated the laminar incompressible and com- 
pressible shear flow along a flat plate by way of series 
solutions in terms of a vorticity number. The same 
problem for the incompressible case has been solved 
by by integral methods. Sakurai" considered 
the case of shear flow along a flat plate with uniform 
suction. Hains!! has presented the formulation of the 
general problem of a laminar boundary layer with non- 
adiabatic rotational free stream suggesting a_ series 
solution in terms of the vorticity number. 

An analogous series solution approach will be used 
here in terms of laminar and turbulent ‘‘vorticity num- 
bers." There has recently appeared to be some con- 
troversy as to the proper boundary and initial condi 
tions to be applied to the basic equations (see reference 
12). In the subject case, the conditions of equilibrium 
of the wake, as suggested by von Karman, make the 
controversial point inessential. 
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Fic. 1. Geometry of the problem 


(1) EguaTION OF MOTION 


The fundamental equations governing the plane 
mixing of an incompressible free stream with adjacent 
still air are here derived for the case in which the free- 
stream vorticity is a constant different from zero. 
Both laminar and turbulent cases will be considered. 


The geometry of the problem is shown in Fig. 1. The 
free stream is described by the equation 
= + (1) 


where m and w are constants. The interaction origi- 
The laminar case will be dealt with 


nates at vy = O. 
first. 


(a) Laminar Mixing 


The fundamental equations are 


UU, + = Uz +t, = O (2) 


The well-known boundary-layer approximation has 
been used. This fact will clearly put some limitations 
upon the relative magnitude of the free-stream vor- 
ticity w. 

The boundary conditions appropriate to Eqs. (2) are 


lim u = + wy, 
(3a) 
lim v = 0, lim 0 
yor 
and the initial condition is 
u(O, v) = uw + wy (3b) 


The second boundary condition, as suggested by von 
Karman, reflects the fact that the wake is not acted 
upon by any resultant force in the y direction. This 
condition makes the controversy about the appro- 
priateness of the first boundary condition’? inessential in 
the present case. 

The first boundary condition can also be written as 


lim 


+2 


u = uo[l + xq] (4) 


where 7 is the Blasius variable 


(5) 


n = V(ito 2vx)' 


and x is the ‘‘vorticity number’ 


xX = (6) 


The width of the wake originated by the interaction 
of a uniform stream mo with still air is given by 


b = R(2vx * k = const. 


Hence, apparently, x is a measure of the importance 
of the free-stream vorticity as compared to the average 
vorticity k mu 6 that would be present in the mixing 
region between a free stream with velocity m and fluid at 
rest. The validity of Eqs. (2) then requires not too large 
values for x. This fact, together with the expression 
lim w given by Eq. (4) suggests a series solu- 
tion of the basic equation in terms of x with coefficients 
that are functions of the Blasius variable. 
Accordingly, it is assumed that’~'* 


y= 


where 7 is the Blasius variable and y is the stream func- 
The velocity components are expressed, to a 


for 


* [go(n) + xgi(n) + +...) 


tion. 
first-order approximation in x, by 


= go + 


= [ngo’ — go) + — 


1/2 
(2 vu) 


Substituting Eqs. (8) into the first of Eq. (2), equat- 
ing terms in the ascending powers of x and neglecting 
terms in (x)* vield the following equations: 


+ gogo’ = O (9) 


— + go + = O (10 


Analogous procedure applied to the boundary conditions 


vields 
lim = 1 lim = 7 | 
na +2 + 
(11) 
lim g’ = 0 
lim lim = | 
+ ro + 
(12) 
lim = 0 | 


Eq. (9), subject to the boundary conditions given by 
Eqs. (11), is the Blasius equation whose solutions are 
already available in the literature.'*'* Eq. (10) is 
linear in g,; and represents, to a first approximation, 
the effects of the outer stream's vorticity. Its solution 
will be given in a later section after having shown that 
the turbulent case can be reduced to the same basic 


equations. 
(b) Turbulent Mixing 


In attendance with the majority of problems involv- 
ing free turbulent mixing, two main hypotheses are 
made here in deriving the pertinent fundamental equa- 
tions. The first hypothesis is that of Reynolds Num- 
ber similarity, which implies that, for sufficiently high 
Reynolds Numbers, the effects of the fluid viscosity are 
limited to the elements of turbulence directly respon- 
sible for dissipation of mechanical energy.” If it is 
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postulated, as in the subject case, that initial boundary- 
layer thicknesses are zero and that the free-stream vor- 
ticity is small compared to the average vorticity pres- 
ent in the wake, no region of the subject flow appears 
to be excluded from the Reynolds Number similarity. 
By this hypothesis, then, the mean viscous stresses will 
be neglected everywhere in the flow field. 

The second hypothesis is that of constancy of the 
eddy kinematic viscosity coefficients over cross sections 
of the flow. This assumption has so far proved to yield, 
in problem of free turbulence, results that fit the ex- 
perimental data with sufficient accuracy. 

By following the lines indicated by Townsend,” it 
can be proved that, even with the acceptance of the 
first hypothesis, no similar functional forms of the 
mean properties and of the Reynolds shear stresses 
are possible in the general case where the free streams 
are of the type “4; = uw» + wy. However if, again, 
some restrictions are accepted as to the order of mag- 
nitude of the vorticity w, a suitable linearization process 
allows the satisfaction of similarity conditions in both 
differential equations and boundary conditions. It 
results that velocity and length scales are the same as 
in the case of uniform free turbulent mixing so that the 
eddy viscosity coefficient is still a linear function of x. 

Thus the basic equation of motion can be written, 
in terms of mean quantities, as 


UU, + = Uy + v, = O (13) 


The boundary conditions are again given by 


lim + wy, 
yo> t+ 
(14) 
lim 


im «= 0, 


v= 0| 


and the initial condition by 


u (0, v) = 


uo + wy 
In terms of a new variable defined as 


the first boundary condition can also be written as 


lim + (16) 


o> +0 
Here ¢ is the ‘turbulent vorticity number” and it is 
given by 


g= (17) 


1/2 
WV EQ Uo 


Series solutions will again be sought in terms of as- 
cending powers of ¢. 

The stream function is expressed, to a first approxi- 
mation in ¢, by 


y = [Go(o) + ¢Gi(c) | (1S) 
The velocity components are accordingly given by 
= Go’ + (19) 


Substitution of these equations into the momentum 
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equation will again show that the functions G, and ¢ 


must satisfy the following equations 
+ GG” = 0 
— Gy'Go! + + = O8 


The boundary conditions are, likewise, 


lim G,' = 1 lim = 
oma ta om + 
2] 
hm Gyo = 0 | 
lin = lim G,; = o? 
(22 
| 


It thus appears that both laminar and turbulent muy. 
ing are governed by the same equations. 

The solution for the linear third-order equation 
vielding, to the first order in ¢ or x, the effects of the 
vorticity, is given in Section (2). 


(2) SOLUTION OF THE Basic EQUATION 
The equation to be solved is 


yer 


+£ 


1 


Zo + 2 


subject to the boundary conditions 


lim =7 lim = 7° 2 
no + 
(24 
hm = 0 


The subject equation is homogeneous and linear; its 
from the 
We must then resort to a numerical 


coefficients are known Blasius solution in 
tabulated form. 
solution. To this end, it proves very convenient to 
determine asymptotic expressions for the solution. 


As » becomes very large, asymptotic values for ¢ 


are substituted into Eq. (23) to get 
+ ng” — gy’ = 0 (25 
The change of variables 
reduces Eq. (25) to 
tH” + (32) — (3 2)H =0 (27 


Solutions of this equation can be given in terms of the 
confluent hypergeometric series F(a, ¢c, ¢). It is known 
that two independent solutions of Eq. (27), analytic 


about ¢ = O and valid for ¢ finite, are given by //; = 
F(a, c, and Hy = Fa —c + 1,2-—69) 
where, in the present case, itisa@ = c = 3 2. It thus 


appears, by recalling the properties of the confluent 
hypergeometric series, that the first solution is simply 
H, = eé; i.e., it is proportional to the inviscid solution 

2 satisfying both boundary conditions at plus 


infinitive. The general solution is then 


(28) 


H = Aj + Ax~'? Fi, (1 2), 


or, through Eqs. (26), 
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dg, df) = Ais 
(29) 
The boundary conditions (24) can be written also as 
follows: 


g(to) =¢, dg, = 1 (30) 


Since the asvmptotie behavior of the confluent hyper- 


eometric series is given by 


F(a, c, T(a)] e (31) 


the asymptotic expression for the solution of Eq. (25) 


valid for n large is 


F[l, (1 2), 


(32) 


As n becomes very large negatively, Eq. (25) can be 
written in the following asymptotic form: 


kg,’ = () 


where = —£(— ©) and is to be considered known. 
The solution of this equation satisfying the second 


boundary condition g;/(— ©) = 0 is 
+ Ay (33) 


With the help of Eqs. (32) and (35) and through any 
standard numerical technique, the solution g; is now 


readily found. 


(3) ANALYSIS OF THE RESULTS 


Insofar as terms in the square of the vorticity num- 
bers have not been considered, the results of the present 
analysis vield the effects of the outer stream vorticity 
on the mixing characteristics to a first approximation 
only. The taking into account more terms in the series 
expansion for x(or ¢) does not in principle introduce 
additional complications. However, such refinement 
of the analysis does not appear justified in view of the 
fact that for higher values of the parameter x the entire 
boundary-laver approach is questionable. 

Since both laminar and turbulent vorticity numbers 
are functions of the x coordinate, the present analysis is 
to be considered valid in the proximity of the origin 
of mixing. The actual range of applicability along the 
x axis is determined, for any given value of the outer 
stream vorticity, by the requirement that the vorticity 
number remains less than one. The vorticity number 
can be considered as ratio between the free-stream vor- 
ticity w and the average vorticity in the wake. As 
this average value decreases with x, a point will be 
reached where the two vorticities are of the same order 
of magnitude and therefore the boundary-layer ap- 
proach breaks down. 

The functions gi, g:’, and [ng:’ — 2g;] are shown in 
Velocity profiles « “,; have been computed for 
several values of the vorticity numbers and are shown 
in Fig. 3. Positive and negative values of the vorticity 
For comparison, the 


Fig. 


numbers are considered there. 
velocity profile corresponding to the uniform stream 
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u = Wl¢e = x = VO) is also shown in Fig. 3. The dis- 
tribution of the nondimensional y component of the 
velocity is shown in Fig. 4 for the same range in the 
values of the vorticity numbers. The effects of vor- 
ticity are readily apparent from these diagrams. 

Besides the velocity profiles in the wake, two other 
characteristic quantities are of practical interest in 
mixing phenomena: (1) the amount of entrained mass, 
and (2) the rate of spreading of the wake. 

The mass entering the wake per unit time and area 
is given by 


pi(v,— ©) = 2x) [(O.S72 + 0.373x] laminar 


pu(x, — ©) = 0.0367 pu0[0.872 + 0.373¢) turbulent 


In deriving the second of the above equations, the 
value of «) has been calculated through comparison 
with the data reported in the Tollmien’s analysis. '® 
The mass flux per unit length, #, drawn out from the 
still fluid by the mixing action up to a station x» is given 


by 
m = 0.S72p(2vuoxy)' {1 + O.214x{, laminar 
m = 0.032puox) }1 + turbulent 
For the case of interaction between the uniform 


stream and still air, the mass entering the wake per 


unit time and area is 


[pv(x, — = 2x) laminar 


[py(vx, — ©)], -9 = 0.032pm% turbulent 


and the mass flux into the wake per unit length is ex- 
pressed as 


= 0.872 laminar 


= 0.032 puexo turbulent 


Thus, at any given station x, the ratio between quanti- 
ties due to mixing of the uniform stream mm and those 
due to the mixing of the nonuniform stream #; are 


v(x, —©)/[v(x, — ~)],-9 = 1 + 0.428x) 
: laminar 
m/([m], = 1 + 0.214x 
v(x, —©)/[v(x, —@©)], -9 = 1 + 0.428¢) 
turbulent 
m/|m), = 1+ 0.2l4¢ 


w 


Concerning these relations, the following remarks can 
be made. 

(1) A positive value of the vorticity in the outer 
stream increases the value of the mass entrained per 
unit time and area as compared to that entrained by 
the corresponding uniform stream mm. The per cent 
increase depends upon the vorticity number and, there- 
fore, for a given vorticity, on the distance from the ori- 
gin of the interaction as well as on the type of interac- 
tions. By considering a value of the vorticity number 
equal to 0.4 as the conservative limit up to which the 
results can be assumed to be valid, the aforementioned 
increase can be as high as 16 per cent. 

(2) Analogous increase is found for the total mass 
flux into the wake up to a given station x». For the 
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Nondimensional y component of the velocity for several 


values of the vorticity numbers. 
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Fic. 6. Application of results of vorticity analysis 
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same upper limit of the vorticity number, the maximum 
increase is of the order of S per cent. 

3) The opposite trend obviously occurs for negative 
values of the vorticity in the free stream. Mass en- 
trainment is less than that ocecurring for the corre- 
sponding uniform stream with velocity uo. 

(4) For a given value of the outer stream vorticity 

these effects are more rapidly and intensively felt 
in the turbulent mixing than in the laminar mixing. 
This behavior stems from the fact that the turbulent 
vorticity number increases linearly with x while the 
laminar vorticity number increases with the squared 
root of x only. Furthermore, in usual conditions, the 
rate of increase is much larger in the turbulent case. 

To make a numerical example, consider a nonuni- 
form outer stream with w mw = 10 ft.-' atx = x = | 
it.and = 150 ft. sec. (wx) v~ 108). In the laminar 
case, the vorticity number variation with xv is given by 


x = 0.014x'"" 


and, in the turbulent case, by 
= 0.367x 


The ratio 7(v, — ©) [v(v, — ©)],, 9 is plotted against 
v in Fig. 5 for both laminar and turbulent cases. It 
can be seen that the effect of the vorticity w is com- 
pletely negligible in the laminar case. The total mass 
entrained into the wake from x 0 tox 10 ft. is 
only | per cent larger than that entrained by the uni- 
form stream with mw = 150 ft. sec. and w = 0. The 
situation is quite different in the case of turbulent 
mixing. At a station x 1.25 ft. the total mass en- 
trained (per unit length) into the wake is already 10 
per cent higher than that entrained by the correspond- 
ing uniform stream at x 1.25 ft. For x greater than 
about 2 ft. the average vorticity in the wake becomes 
so low that the present analysis is very likely to be 


defective. 

It would be interesting to determine the rate of 
spread of the wake in the subject cases. It is well 
known that the arbitrariness in the definition of a 
wake width is usually overcome by taking its outer 
edges to be the loci of points where the velocity at- 
tains a given per cent value of the outer streams. 
When the outer streams are uniform the flow field is 
sunilar and, consequently, the wake width thus defined 
is a well-determined function of x. When a constant 
vorticity is present in the outer stream the velocity pro- 
files are obviously no longer similar. A well-defined 
variation of the wake width with x cannot be deter- 
mined. It can, however, be said that the wake width 
will certainly be no longer a parabolic function of x 
(laminar case) or a linear function of x (turbulent case). 
Other interesting remarks can be made by considering 
the function g;’._ From Fig. 2 it appears that g;’ is neg- 
ative for (—n) sufficiently large. This obviously im- 
plies that for a positive value of the vorticity number 
the vanishing value of the velocity is attained faster. 
On the other hand it appears that, since g:’’ approaches 
the value 1 from below, the function g,’ approaches its 
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asymptotic value from above. That is, for 7 large, 
gi >. Thus, in comparing the wake width due to the 
stream 1) uy + wy with that due to the stream mM, it 
can be said that, at least initially, at any given station 
x the Feight of the dissipative region is smaller for 
w > O (and, conversely, larger for w < O) than that 
This result cannot be extrapolated to large 

Indeed, for x (or ¢) greater than about 


forw = 0. 
values of x. 
0.4, an overshooting in the velocity profiles is noticed, 
which, very likely, indicates the failure of the approach. 
A similar breakdown of the analysis is to be expected 
for large negative values of the vorticity numbers. 

All the above results can be summarized by saying 
that, in practical situations, any effect of the constant 
vorticity present in the outer stream is negligible over 
a considerable distance from the origin of the interaction 
in the laminar case. This effect is instead rather read- 
ily felt in the turbulent case due to the rapid rate of 
increase of the turbulent vorticity number with x. 

The present results are of interest from another point 
of view also. As there are practical cases in which the 
properties of the free streams are given, so there are 
cases when, instead, certain mixing behaviors are pre- 
scribed and or desirable. For instance, in injectors 
and augmentors, which involve the mixing of two 
streams of different momentum, it is important, from 
a practical point of view, to achieve large mixing rates 
so that the device can be made as short as possible. 
On the other hand, in the use of a low energy stream 
interposed between a surface to be protected and a 
high energy outside stream, it is desired to reduce the 
mixing between the two streams so that the length of 
surface protected for a given mass flow of the low 
energy stream is comparatively large. Finally, prob- 
lems connected with combustion, ignition, and flame 
stabilization furnish other evident instances of flow 
fields where the controlling of the rate of exchange of 
properties between two streams is essential for the efli- 
ciency of the process itself. Hence, the question arises 
as to whether and how the presence of vorticity in the 
external stream could be used to achieve these ends. 
An answer comes immediately from the present analysis 
when applied to streams having constant vorticity 
over a finite y length (which is obviously assumed to 
be much larger than the dissipative region's height). 
Consider the flows shown in Fig. 6. <A stream of ve- 
locity u.. is to mix with a quiescent region Fig. (a). 
If it is desired to increase the rate of mixing, a stream 
of higher velocity fluid could be injected upstream of 
the origin of mixing so as to obtain a velocity distri- 
bution, as shown in Fig. 6(b). Indeed the ratio of 
mass entrained up to a station x = » for the flow con- 
figurations (b) and (a) would be 


= (uo/u.)/? {1 + O.214x}, =, 


(m)»y (M)a laminar 


= {1 + 0.2l4¢;,-,, turbulent 


(M1), (Ma 


and thus would be greater than one if % > u., since, 
in the worst case, the term in braces would amount to 
Analogously, if the mixing is to be retarded, 


0.92. 


{ 
| 
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then a screen or vorticity generator can be located so 
that, at the origin of mixing, the velocity profile is as 
shown in Fig. 6(c). It is of interest to point out that, 
in the case where a low energy stream is injected over 
the surface of a hypersonic wing in order to insulate it 
from the high energy flow, it is simple to realize the de- 
sired vorticity distribution for reducing the mixing. 

The above discussion is indicative of the practicabil- 
ity of mixing studies involving vorticity in the free 
streams and suggests future investigations as to the 
effects of vorticity on temperature profiles for low-speed 
flow and the extension of the analysis to more general 
cases of incompressible and compressible mixing. 


CONCLUSION 


The incompressible laminar and turbulent plane 
mixing of a nonuniform constant vorticity stream with 
fluid at rest has been solved to the first approximation 
in the vorticity numbers. It has been found that, 
when the characteristics of the mixing due to the non- 
uniform stream “4, = Wp + wy are compared with those 
due to the uniform stream w, the effects of a vorticity 
in the outer stream are mostly negligible if the inter- 
action is laminar, whereas they are sizable if the inter- 
action is turbulent. In this latter case, for instance, 
the increase in total mass flux entrained up to a given 
station x may become several per cent larger than the 
one drawn out by the corresponding uniform stream. 

The different effects of vorticity on laminar and 
turbulent flows are due to the different variation of 
vorticity numbers with «. Laminar vorticity number 
increases with the square root of x and turbulent vor- 
ticity number with «the rate of increase being, in 
usual instances, much larger for the latter case. 

Application of the results to streams having constant 
vorticity over a finite y length has shown that the 
presence of vorticity in the outer streams can consti- 
tute a valid means to control the rate of intermixing 
of two streams. This conclusion warrants extension 
of the analysis to more general cases in view of the 
variety of problems in which the rapid transferring of 
some properties from one stream to another is of sub- 
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stantial importance for the efficiency of the process 
themselves. 
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Special Summer Course in Heat Transfer 

Heat Transfer Laboratory, University of Minnesota | 
if 
The Heat Transfer Laboratory, of the Mechanical Engineering Department, University of Minnesota, announces a 
special course in Heat Transfer to be given from August 7-15, 1958. The program will emphasize modern developments 

in heat transfer and is aimed at engineers and scientists working in the aircraft or missile field or in the nuclear power field. 
Such topics as Mass Transfer Cooling, Heat Transfer at High Altitudes, Radiation Properties of Gases at High Tempera- 


Hartnett, Mechanical Engineering, University of Minnesota. 


tures, Reactor Heat Transfer, Two Phase Flow, and Boiling will be covered. In addition to the staff associated with the 
Heat Transfer Laboratory operating under the direction of Dr. E. R. G. Eckert, guest lecturers will participate. 
Tuition for the course will be $180 and the enrollment is limited. Jnformation may be obtained by writing to Dr. James P 
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On Vibrations of Conical Shells 


G. HERRMANN* ann I. MIRSKY** 
Air Research and Development Command, USAF 


SUMMARY 


Frequencies of truncated conical shells are determined using 
Rayleigh-Ritz procedure. Sinusoidal mode shapes are as- 
sumed, which satisfy the equation of motion of a corresponding 
evlindrical shell. The approximate values should therefore have 
90 id accuracy for small semivertical angles, or small ‘‘conicity.’ 

It was found that for short shells the conicity lowers the fre- 
quency somewhat, while for long shells the frequency increases 
appreciably with conicity This influence is strongest if the 
number of circumferential waves is 38. As may be expected, the 
frequency of thin shells and the lowest frequency are more sensi- 


tive to changes in shape. 
It is shown further that the equations of motion may be inte- 
grated in terms of known tabulated functions in several special 


cases. 


INTRODUCTION 


—Pyngeeenges ATTENTION Was given in recent years 
to stress and stability analysis of conical shell 
elements (see, for example, reference 1 and the bibli- 
ography to reference |) while the dynamic behavior of 
such shells remained relatively unexplored. The most 
recent paper that has come to authors’ attention on 
this problem was the study by Federhofer.* In this 
investigation, Federhofer determined the natural fre- 
quencies of complete and truncated conical shells of 
arbitrary semivertical angle and with several different 
boundary conditions. The frequencies were determined 
by application of a Rayleigh-Ritz procedure, assuming 
the displacement components to be polynomial fune- 
tions of the coordinate along the axis (or the generator) 
of the shell. Since it is well known that some mode 
shapes of cylindrical shells are trigonometric functions 
of the coordinate along the axis (see, for example, refer- 
ence 5), Federhofer’s results cannot be expected to yield 
accurate values for the corresponding frequencies in 
the limit for cylindrical shells. 

In contrast to the aims of reference 2, the purpose of 
the present study is to establish the influence of slightly 
conical shapes on the frequency of free vibration, as 
compared to frequencies of corresponding cylindrical 
shells of same thickness and length. Therefore, in 
applying the Rayleigh-Ritz procedure, the same dis- 
placement functions will be used as those which are 
known to satisty the equations of motion of cylindrical 
shells, and only small semivertical angles (small 
“conicity’’) will be considered. 

The first section contains a brief discussion of the 
assumed deformations, and in the second section the 

Presented at the Missile Structures Session, Twenty-Sixth 
Annual Meeting, IAS, New York, January 27-30, 1958. 

\ssociate Professor of Civil Engineering, Columbia Uni- 
versity 

* Development Engineer, Foster Wheeler Corporation. 


energy expressions are derived. The third section deals 
with a general form and several special forms of the fre- 
quency equation, obtained by application of a Ray- 
leigh-Ritz procedure. A discussion of numerical re- 
sults is presented in the fourth section. 


tion contains a discussion of certain special solutions of 


The last sec- 


the equations of motion, for axial symmetry, and ne- 
glecting bending effects. 


DEFORMATIONS 


We consider a conical shell of median radius R, 
length / and semivertical angle a. 
of the shell are referred to an x, 6, 2 coordinate system 
x is the ‘‘generator coordinate’’ measured 


The mass elements 


(see Fig. 1). 
along the generator, starting at the midpoint .1/ in the 
direction away from the vertex. @ is the circumferential 
coordinate measured in planes perpendicular to the 
These two coordinates locate points 


zis the “normal 


axis of the cone. 
on the middle surface of the shell. 
measured from the middle surface in the 
The components of displacement 


coordinate,” 
outward direction. 
in this coordinate system will be designated by w,, 
u,, and “,, respectively. 

To find the components of strain in this coordinate 
system, we observe that the square of a length element 


ds is 
ds? = dx? + (R+ xsin a + 2cos a)? dé? + (1) 
This expression is readily obtained from the corre- 
sponding relation in cylindrical coordinates 7, 6, 7 = 
R + 3, 
ds? = + (R + 2)? + d® (2) 

using the transformations 

= xcosa — a) 

6=80 (3) 

= x sin a + 2 cos 
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The components of strain in a general, orthogonal, 
curvilinear coordinate system are expressible in terms 
of coefficients in the expression of the element of length. 
Using the relations, as given for example by Love,‘ we 
obtain the following components of strain: 


= On./Ox 

Cog = [(Ou, + sin a + cos a| + 
(R+xsin a + 2 cos a) 

= On,/OZ 

Yo: = (Ou,/0z) + 06) — u, cos a] + (4) 


(R +x sin a + 2 cos a)t 

(Ou, Of) + (Ou;/Ox) 
Yro = (Ou, Ox) + }[(Ou,/ 0) — Uy sin a] + 
(R+xsin a + z2cos a)} 


| 


It will be assumed, similar to previous work on shells,° 
that the displacements “4, are approximated by 
displacements @, respectively, whose dependence 
on the z coordinate is specified in the form 


= (Ou/Ox) — 2(0°w/dx?) 
Pog = (R + xv sin a)] + 
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( = u(x, ¢,t) + 2 
(x, = w(x, 9, | 


where (is the time. 


It is noted that the functions , v, w represent dis 


placements of the elements of the middle surface, whij 
y, and y, represent angles of rotation. 
the cross sections remain perpendicular to the deform, 


middle surface—that is, neglecting the effects of trans 


verse shear deformation— we put 


As a consequence, ¥, and y, are expressible in term 
of the displacements of the middle surface as 
y, = —Ow/Ox | 
= [vcos a — (Ow O¢)| (R + x sin a)f 


The expressions for the components of strain unde 


these assumptions are 


sin a(Ow (R + x sin a + cos — 


+ x sin a) (R + x sin a + 2 cos a) 


oo 
= Voz = ¢:, = O 
= [(Ov Ov) (R+ x sina +scos a) (R+ sin a)] — 


(R+ xsin a) (R+xsin a + cos a)] 


(2R + 2x sina + scosa) + 


[v sin + xsin a+ scos a) (R +x sin a)*] + 


[z sin a (Ow/O¢) (2R + 2x sin a + scos a) (R+ x sin a)? (R + x sin a + 2cos a)! + 


ENERGY EXPRESSIONS 


The strain energy of a volume |” 


WW = (1/2) (9) 


can be expressed in terms of components of strain using 
the stress-strain relations 


@ [E v?) | + veges) 
Foe [E v*)| — very) (10) 


ll 


Ore 


which, in turn, can be derived from Hooke's law of 
three-dimensional isotropic elasticity by eliminating 
e,, and neglecting ¢,,.. The application of this procedure 
in thin shell theory was justified in reference 6. FE 
designates Young's modulus and y Poisson's ratio. 
Noting that the element of volume in the present case is 


(R + x sin a + 2cos a)dxd odz 


the strain energy of the shell of thickness /: is 


cos a@ h/2 
W = (1 » [E/(1 — X 
—l 2cosadO0 —h/2 


+ + + [C1 — X 
(R+ x sina + 2cos a)dxdg¢dz (11) 


The kinetic energy 7 of the shell, in case rotatory in- 
ertia is neglected, is given by 


[(Ou Od) (R + x sin a + COs 


cos aJ —h/2 


(R+xsin a + scos a) dxdddz (12 


where p is the mass density. 


FREQUENCY EQUATION 


The familiar Rayleigh-Ritz method is used for the 
determination of the frequency equation. In applyim 
this procedure, the displacement components are as 
sumed to be of the form 


u = L’ sin (mx cos a /) cos noe™™ 
u = V cos (rx cos a /) sin noe > (13 
w = W cos (rx cos a 1) cos noe" 


where n is a positive integer indicating the number oi 
circumferential waves. 
symmetric deformations. It is noted that these dis 
placements satisfy the boundary conditions 


ou/Ox =v=w=O0 at x = +/ 2cosa (It 


It is expected that these displacements will yield rea 
sonably accurate frequencies for slightly conical shells, 
say 0° < a < 15°, since they satisfy the equations 0! 
cylindrical shells, a = 0. 


Assuming tha: 


n = 0 corresponds to axially } 
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Substituting the forms (15) into expressions (11) and 2a, ay by 
12) with the aid of the strain-displacement relations als 2), Go| = © (15) 
and minimizing the total potential + 7 with be C2 2c 


re 


spe 


ct to the amplitudes LU’, I’, and IV, the frequency 


It where the following notation has been adopted: 
ation, in determinant form, results as 
juation, 1 ’ 


l R = @ = 
h/R=m w,? = E/pR*1 — v? 


a, = cos 2j) + [sin? + (1 — 2] (0, + m6; cos? 12j7) — (@j 2 cos a) 
= + sin? a(l — v) 2] — 2cos a) + [((1 — v)m? cos? X 
sin a cos a + sin? + cos* + cos a) 
c = (Am? 1272) (n? — cos? a)? + (rtm? cos* a 247°) — (a@j/2 cos a) + 6:(1 + 67") cos? a + 
— v)n*@ym? sin? a@ 677] + 127°) [sin? + 2(1 — v)n*\cos* + 
sin a cos a[(3 — 2v)n? — cos? a) 
»=al(l + + (38 — sin a] 2 
= vr cos a + cos* a 127?) + 6; sin cos + a — (3 — v)n® 2] sin a cos a 1277; — 
\@sm? alsin? a + (1 — v)n® 2) cos? 
= a(l — v)n 47?] + cos a(1 + a@ 127?) + sin? a cos a(l — v)n 477] + 
sin a cos? a(l — v)n 
= [1 (2 sin a)] jlog,[(2 + jtan a) (2 — j tan a)] — [C24 cot aj + — cot a j — X 
cos (27 cot a j) — cot a j + — cot a j — sin (27 cot a Dt 
= [1 (2 sin jlog,[(2 + j tan a) (2 — j tan a)] + [C24 cot a j + — cot a j — X 
cos (27 cot a j) + [S\(27 cot aj + — cot a j — sin (27 cot a 
= (1 sin a) } [S24 cot a j + — cot a j — cos (27 cot a j) — cot a j + 
cot a j — sin (27 cot a 
4, = cot? a} log,[(2 + j tan a) (2 — j tan a)| sin a — 
= —0, + seca j[(1 7?) — (tan? a 
= cot? a — sin a[(1 j*) — (tan? 4)] 
6. = —26,tana 
A, = 70; cot a + sec a [(1 7?) — (tan? a 4)} (16) 
The above 6's contain the functions C,(x) and S;(x) 
which are the cosine and sine integrals defined as The frequency equation of the membrane theory is 
a obtained by setting m = 0 in expressions (16). For 
C(x) = -| (cos ¢ t) dt non-axially symmetric motions it is still given by 
(17) Eq. (15); however, the constants a), ete., are now 
Si(x) = (sin ¢ dt 
0 a, = cos 2j) + — 24+ 
and stem from the integration of the total potential sin? a] — (aj 2 cos a) 
over the length of the shell from = —/ 2 cos a@ to b, = — ») sin* a2 + 
1/2 cos a. (1 — pv) (r* cos a 47) — (aj 2 cos @) (21) 
= 62 cos? a — (aj 2 cos a) 
SPECIAL CASES dy = [n(1 + v)r 2! + n(3 — v) (6; sin a 2 
by = vr cos a + sin cos a 
Several special cases may be obtained directly from iu te ee 


the determinant (15). For the case of axially sym- 


metric motions the frequency equation (15) simplifies to 


As before, the frequency equation for axially sym- 
metric motions is 


2a, by 
“ 0 (18) 
by 26 40 = 0 
where dj, 62, ¢ are obtained by setting n = 0 in the ex- where the constants a, ¢:, and }, are obtained by set- 
pressions for d;, b:, (4, respectively. It governs the two ting = 0 in the set of expressions (21). 
flexural modes. The torsional mode is obtained by The torsional mode is given by 
setting 
bh = 0 (22) 
IW) oV = 0 (19) 
or 
lhis yields the result 
— v) 2] + [m? cos? a(l — v) j?] X (1 — v)@ sin? a + (1 — v) cos a — 
(aj cosa) = 0 (23) 


[0-7 sin @ cos a + 4 sin? a + 672 cos? a + 


(2277 m? cos a)| = aj 2cosa (20) Finally, it is of interest to perform the limiting process 


! 
dz (12 
5 
(13 | 
nber 0 
axially | 
‘se. dis- | 
a (14) i 
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a — 0 and obtain the frequency equation of a cylindrical 
shell as a special case. As a — O, the functions C,(x) 
and §;(x) appearing in expressions (16) have to be ex- 
panded in asymptotic form. The asymptotic expan- 
sions emploved are 
~ 2) — P(x) cos — Q(x) sin x| (24) 
~ P(x) sin x — Q(x) cos x f 
where 
P(x) ~ (1/x) (21/x*) + (4!/x*) 
Ole) ~ ~ +... 


With these asymptotic forms the constant 6’s can read- 
ilv be evaluated. The simplified expressions become 


A = & = 2, 0; = Ol (25) 
and the constants etc., reduce to 
ay = (9? + (1 [1 + (m? 12) 2} — aj?) 
by = + (1 — v) [1 + (m? 4)] 2} —aj?) 
= (rim? 247°) + (m2n?x? + 
+ m? (n? — 192/12] — (aj 2)} 2 (26) 
a, = (1 + 
bo = ve + 127? — (1 — /24 
= + — v) 24] 


The frequency equation (15, using expressions (26), 
may be verified to be identical to the one based on the 
Fliigge equations’ of the bending theory of cylindrical 
shells. 


NUMERICAL EVALUATION AND DISCUSSION OF RESULTS 


The numerical evaluation of the frequency equation 
(15), based on the bending theory of conical shells, was 
carried out for three values of the semivertical angle 
a—namely, a = 5°, 10°,and 15°. This angle a can be 
taken as a measure of “‘conicity’’ of the shell. Further, 
the number of circumferential waves ” was taken from 
0 to 6 and the ratio of the length of the shell / to the 
mean radius R from 1 to 10. Two thickness to mean 
radius ratios were considered—namely, m = 1/30 and 
m = 1/100. 

Since the main purpose of the present investigation 
is to establish the influence of conicity, it appeared 
most useful to represent the frequency of a conical 
shell (a # O) as a ratio to the frequency of the corre- 
sponding cylindrical shell (a = 0), rather than as a 
ratio to some fixed reference frequency. 

Frequencies of thin cylindrical shells were tabulated 
in a very complete fashion by Baron and Bleich.* 
In these tables the frequency w was represented in the 
form 


w = K RV(G/p) [1 + m(K?2/K?)] 27) 


and the factors K and K? K? were tabulated for differ- 
ent m and for different ratios 1’R. It was shown that 
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NON: TORSIONAL / 


je 


Fic. .2. 


the second term under the square root contributes only 


to the lowest frequency for n > 2. 

To determine the ratio 2 of the frequency of the cone 
to the frequency of a corresponding cvlinder, we have 
to multiply Wa = w w by 


= VE pR*(1 — 
and divide by 
K RV(G p) (1 + m?) (KR? K? 


The result is, for vy = 0.3, 


[1.69 KV 1 + (@ (28 


Fig. 2 shows the plot of 2 versus LR for the lowest 
nontorsional mode, ” = 0, for three values of conicity 
namely, a = 5°, 10°, and 15°. The abscissa itself is, 
of course, 2 for a = 0°. The frequencies are independ- 
ent of m. It is seen that for short shells the 
conicity lowers the frequency, while for long shells the 
frequency is increased and is strongly dependent on a. 
For the torsional mode the plot of Fig. 3 is essentially 
the same, except that numerically the influence of 
conicity is weaker. For the higher nontorsional mode 
n = O, it is still weaker, as seen from Fig. 4. 

The frequency of the lowest mode for 7 = | is plotted 
in Fig. 5, and again the same qualitative dependence on 
L Rand on ais noted. 

The frequency of the lowest mode for 7 = 2 becomes 
dependent on the thickness to mean radius ratio R. 
This dependence is illustrated in Fig. 6. As may be 
expected, a thinner shell is influenced more strongly 
by conicity than a thicker shell. 

To illustrate, for the lowest mode, the influence of 
the number 7 on the frequency as a function of a, plots 
of Fig. 7 have been prepared. L R was chosen as 
7 and m as 1/100. It appears that the influence of 
conicity is strongest for 2 = 3. 

Finally, Fig. 8 indicates the influence of conicity on 
the three modes, as a function of LR, for » = 3 and 
a = 10°. As may be expected, the influence is much 
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stronger for the lowest mode and weakest for the high- 
est mode. 


DERIVATION AND ANALYSIS OF EQUATIONS OF MOTION 


In addition to determining approximately the nat- 
ural frequencies of a conical shell, an attempt was 
made to solve the governing equations of motion rigor- 
ously. However, even in the simplest case of axially 
symmetric motions and using membrane theory, these 
equations of motion could be solved only under certain 
supplementary assumptions, as will be discussed be- 
Therefore, no numerical calculations were car- 
ried out. The equations of motion themselves will be 
briefly derived first, using Hamilton's principle. 

In case of axial symmetry, the strain-displacement 
relations (4) reduce to 


= 


low. 


Cog = (U Sin a + woos a) + (29) 
(R + scos a+ sin a) 

where uw and w are functions of « and ¢ only. The 

strain energy, per unit of volume, reduces to 

21V* = o,,(0u OX) + sin a + wecos a) + 
(R+scosa+xsina)] (30) 


The total strain energy of the shell is obtained by in- 
tegration of expression (30) over the cross-sectional area, 
the element of area being 


(R+scosa+xsin 


and over the length of the shell, from —/ 2 to / 2, with 
the result 


rf [.V..(R + x sin (Ou Ox) + 


= 
Nos (Wcosa + usina)|dx (31) 
The stress resultants .V,, and .V,, are defined as 
cos a 
Nia Grz{1 + - ) dz 
—h/2 R+xsina 
(32) 


h/2 | 
Noo = j | 
—h/2 


The variation of the strain energy expression is given by 


= 


: (} —(0/0x) (R + x sin + 
Veg sin a} 6u + Nog cos a bw) dx + 
(x sin a + R) bw]! » (33) 


The total kinetic energy of the shell is 


lr h/2 
T = (1 » f + w?) X 
-l/2/70 —h/2 


(R + zscosa + xsin a) dxd¢dz (34) 


and its variation, neglecting rotatory contributions, is 


1/2 
= (1, 2) f ph(R + x sin + w*) dx (35) 


-1/2 
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Hamilton's principle, 


bs 
sf = (36 


is now applied to deduce the stress equations of mo. 
tion, together with the appropriate boundary condi. 
tions. JL is the Lagrangian 7° — IT’. 
put into the form 


le 


which becomes, after substitution of Eqs. (33) and (35 


Eq. (36) may be 


and partial integration with respect to time, 


—ph(R + x sin — Ny, sin a + 
1/2 


(0 Ox) [.V,,(R + x sin a)]jéu + 


[—ph (R + x sin a)w — cos aliw) dtdx — 


fe 
(R + x sin + 
ty 


[ph(R + x sin + wéw) dx (88 
—1/2 


This expression must vanish identically for arbitrary 
variations 6” and éw. Hence, the stress equations of 
motion are obtained by setting equal to zero the factors 
of 6u and éw in the first integrand of Eq. (3S), with the 


result 
(0 Ox) }.N,,(R + x sin a)} — N,, sin a = 

ph (R + x sin a) 
ph(R + x sin a)w 


(39 


— Ny COS a = 


The last two lines of Eq. (38) could be shown to contain 
the initial and boundary conditions, which, however, 
will not be discussed at present. With the aid of stress- 


strain relations 


Orr = [E/(1 — + vege) t 


40) 


we derive the appropriate stress resultant-displacement 
relations, using Eqs. (29) and (32). They are 


N,, = E,(Ou Ox) + E,v[(u sin a + w cos a) + 
(R + x sin a) (41 
= FE,[(u sin a + weos a) + 
(R + x sin a)| + E,v(Ou Ox) 
where E, = Eh/(1 — v?). 


The displacement equations of motion then take the 
form 


E, 0) [(R + x sin a) (Ou Ox) + 
v(u sin a + Ox} — 
E,[(u sin a + wcos a) sin a + 
(R + x sin a)| — FE, v sin a(Ou Ox) = 
ph(R + x sin a)ti 


(42 


— E,|(usina+ weos a) cosa (R+xsin a) |— 
E, v cos a(Ou Ox) =ph(R + x sin aw 


These equations are identical to those presented in 


W 
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dt = 0 (37 
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reference 2. Assuming periodic motions— that is, 
u(x, t) = U(x)e™ | 
ict (43) 
w(x, t) = 


the equations of motion reduce to ordinary differential 

equations 

+ Ul’ sina + 

+ x sin a) E,| — 

(sin? a (R + x sin a) U + rcos al” — 
lsin acosa W (R + xsin a)] = 0 


R+xsin 


(44) 
‘+ [sin acosal’ (R + xsin a)| — 
[phw*(R + x sin a) E,] — 
[cos? a (R + x sin a)}} = 0 


y COS @ l j 


where primes denote differentiation with respect to x 

A solution of this set of equations in terms of tabu- 
lated (or previously defined) functions was not found, 
except in a number of special cases in which additional 
simplifying assumpuions have to be introduced. A 
brief discussion of these cases follows. 


Case (1) p Proportional to (R + x sin a)? 


Using this assumption and performing the transfor- 
mation 
y=R+xsinea (45) 
a single equation in l’ is readily obtainable, the result 
being 
+ av? cos? a) + 
v(1 + av? cos? a) (dU dy) + 


(6—acos?a)l = 0 (46) 


where 


a = E, (khw* — E, cos* a) | 


b = (khw? — E, sin? a) E, sin? af 


and & is the factor of proportionality 
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p = k/y? (4S) 


The differential equation (46) is homogeneous and its 
solution is of the form 

U = + By” (49) 
where A and B are arbitrary constants of integration 
and is 


8? = (b — a cos* a) (1 + a) (90) 


In terms of the coordinate x the displacements l° and 
IV take on the final form 


l’ = A, cos [8 log (R + x sin a)| + 
A, sin [8 log (R + x sin a) | 
W = $(A, + vBA2) cos (91) 


[8 log (R + x sin | + — sin 

[8 log (R + x sin asin a cos a 
The frequency equation results from the satisfaction of 
the boundary conditions 
U=0 at x = +/ 2cosa (52) 

vielding 
1+ (/tan a) 2R}j = 

ie, 1,2... Goa 


8 log }[1 — (/ tan a) 2R} 


Case (2) Transverse Inertia Neglected 


With this assumption the resulting differential equa- 
tion in L’ is given, again using transformation (45), by 


+ (dU dy) + 


sin? a(1 — =O (54) 


where n> = w*ph E, (55) 
Eq. (54) is a Bessel equation of zero order, and its solu- 
tion is 
/ 
U = AJj(ny/sin aV 1 — v*) + 
BY\(ny sin aV 1 — 


(56) 


In terms of the original coordinate x, the mode shapes 
are 


U = AJy[n(R + x sin a) ‘sin 1 + + + x sin a) sin aV 1 — v?] 


= (nv tan sinaV1 — {AJ [n(R + x sin ‘sin aV1 — (57) 
BY,|[n(R + x sin a) sin aV 1 — — tan a} AJo[n(R + x sin a) sin 1 — + 
BY,[n(R + x sin a) sin aV1 — 
The frequency equation is obtained by setting either L’ or .V,, equal to zero at x = +/ 2 cos a. 
Case (3) x sina/R<1 


This restriction applies to short shells and to shells with slight conicity. 


Lettings = x sin a Rand omitting terms 


with squares and higher powers of z, the equation on L” becomes 


(a + bs)U" + (c + dz)U’ + (e+ fz)U = 0 


(58) 


where primes denote differentiation with respect to s and the constants are designated by 


a = (n*R* — cos? a) [n?R? — (1 — cos? al 
b = 2(n*R* — cos? a) [3n?R? — (1 — v?) cos? a] + 2n?R2y? cos? a 
= a — 2n*R*y? cos? a 


| 
ix (38 
39 


By a further transformation, Eq. (5S) can be reduced to 
a confluent hypergeometric differential equation (see 
p. 434 of reference 8). 
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to obtain a theoretical check of the experimental critical 
times of 2.00 and 1.80 hours. 

A rather comprehensive series of compressive creep 
tests on V-groove edge supported 2024-T3 aluminum 
alloy plates is presented in reference 7. The results of 
this investigation clearly indicate that there may be a 
significant difference between buckling and crippling of 
flat plates under certain creep conditions. 

The plates tested had 6h values of 20, 30, 45, and 
60. At temperatures of 400°F., 450°F., and 500°F., 
the specimens with } / of 45 and 60 buckled upon the 
initial application of load and required several hours of 
creep before failure occurred. On the other hand, the 
b h = 20 specimens experienced creep crippling at 
strains corresponding rather closely to the theoretical 
critical strain computed from Eq. (6) at the three tem- 
perature levels used. The theoretical critical strain for 
b h = 30 specimens was quite conservative as compared 
to the experimental strains at creep failure. 

From this brief review, it appears that the critical 
strain approach may be of value in predicting creep 
crippling of flat plates of relatively small b/h values. 
The present approach is of no value for predicting the 
creep crippling of flat plates that buckle upon the initial 
application of load. This is essentially a different prob- 
lem and is treated in reference 8. 


(7) CONCLUDING REMARKS 


The critical strain approach presented in reference | 
apparently meets with some initial success in predicting 
the critical times of 3003-0 aluminum cylinders under 

zs compression and torsion. Before the approach can be 
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= 4n*R? — (1 + cos? a] + — cos? a) [n?R? — (1 — v?) cos? a] 
[(m?R? — sin? a) (n*R* — cos* a)? — 2n?yR? sin? a cos? a — sin? a cos? a(n?R? — cos? a)] sin? a 
[4n°R? (n?R? — cos? a) (n?R? — sin? a) + 2n?R? sin? a cos? a + 
2n*R? (n*R* — cos? a)? — 4n?R? sin? a cos? a] sin? (56 


* Love, A. E. H., A Treatise on the Mathematical Theor, 
Elasticity, 4th Edition, Dover Publications, New York, 1944. 
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A Critical Strain Approach to Creep Buckling of Plates and Shells 


(Continued from page 434) 


accepted with confidence, a considerable experimental 


investigation requiring other materials with a wide | 


variation in stress and temperature levels is required. 

The approach does not apply to flat plates that buckle 
upon the initial load application and carry loads in the 
post-buckled state up to the time of creep crippling 
The critical strain approach, however, does appear to 
have predictive value for plates of b h < 20 according to 
the small amount of test data available. 
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Flows in Partly Dissociated Gases* 


Mantred Heil 

Institut fur Theoretische Physik der Freien Universitat Berlin, 
Berlin, Germany 

February 6, 1958 


| gman the laminar boundary-layer equations for a partly 
dissociated gas have been solved for a stagnation point flow. 
The object of these studies was, particularly, the dissociation non- 
equilibrium. Here we want to give equations specified in refer- 
ence 1 for a rather simple case of a partly dissociated diatomic 
gas Apo, in which the influence of diffusion and the transition to 
the corresponding equations for the equilibrium is easily recog- 
nizable 
If we introduce the degree of dissociation 

a = mn,/p, p = mn, + mons (1) 
and neglect the thermodiffusion flow, the following stationary 
boundary-layer equations result for axial flow round a symmetric 
(infinitely long) body of revolution, where we have used the nota- 
tion of reference 2: 


(prvz).2 + (prv,).y = 0 (2 


+ plyay = (1 + ayy} oy + 
(pken/2) {1 — — [4Pa?/K,(1 + @)}} (8 
+ = —P.z + Py = 0 (4 

+ — = + (ATiy),y + 
| (1 + ry (5) 


Dy» = nDyp, w, = (5/2kT + Ep — P = nkT (6) 


* The research discussed in this paper was sponsored by The Martin 


Company. 
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= (1/m,) [((5/4) (1 + + (1/2) (1 + (aEp/2)| 

(7) 
One can calculate values for the velocity of reaction ken |where 
we set approximately k; = ke» = ke (compare reference 2)], for 
example, from the reaction cross sections given in reference 6. 
To the parabolic system |Eqs. (2) to (5)| enter still boundary 
and initial conditions for the unambiguous determination of 7,, 
v,, 2, a. Application of the similarity transformation given by 


Lees’ and consideration of the relations valid for stagnation- 
point flow leads then to a system of ordinary differential equa- 
tions.'. We can write the energy equation (4) in still another 
form by means of 
= Cor Qyr (8) 

From Eq. (7) it follows then 

Cra = (5/2) (R/m)a + (1 — a@)Cy, Cor = wp/me (9) 
If we further introduce 


= UCpa/dA, = wl + = = 
Schmidt Number (10) 


Le = = = Lewis Number (10’) 
the energy equation can be written 


+ plyly — = 
+ [(u/cadiy + — 1)],, (5’) 


and the diffusion equation (3) becomes 


ptra,s + ptya,y = [(u/op)a,y|,y + (pken/2) X 
—a — [4Pa? K,(1 + a)|! (3’) 


Now we want to go over to equilibrium. For kgm ~ © we get, 
since 
nh, = (pkpn 2) = |4Pa? K,(1 + a)|} < @ 


tpa,?/(1 —a,*?) = K,( T), a= al P, T), 


and Eq. (5’) results in 


puri,s + PUylyy + 
+ Cprayr(Le — 


Now we have 
= CyT,, + Cp = Cra + Coranr (12) 


Further, for the equilibrium constant we have 


(d log K,)/dT = w,/kT? (13) 
Together with Eq. (11), it follows that 
ayr = [afl — a,?)/2kT?|w, (14) 
Therefore, 
= Cop = a1 — a,?)w,?/2kT?*me (15) 


and, thus, 
+ Le — 1)T,, = 
(A + =r = (uC,/o)i,, (16) 
Ap = = = Prandtl Number for the 
dissociated gas A» in dissociation equilibrium (17) 
For this we can write also the Schmidt Number op = uwCpp/Xo. 
The energy equation becomes 
+ PUylyy —0:P,2 = BU [(u/o sy (5") 
The diffusion equation (3’) degenerates into a defining equation 
of the reaction rate mr; and Le becomes superfluous. 
The relation 
Ap = [pDiza,(1 — ay?)/2mkT|w,? = 
[PDi2a4(1 — ag*)/2R?T*] (w,/me2)? (18) 
had been already substantially deduced by Nernst’ (in 1904!) 
(up to a factor 2/(1 + a,)*, which results from a somewhat dif- 
ferent assumption for the diffusion flow). It is really unintelli- 


gible that the simple relations C, = Cpa + Cppand X = \ 4 
for a dissociated diatomic gas 12 in dissociation equilibrium hg 
not been used for a long time in the boundary-layer theory 3 
although they are given in reference 8 and easily follow fr 
reference 3. 

If we consider also the thermodiffusion aloug with the diffusj 
which comes from the differences of concentration, we get, j 
stead of Eq. (18), 


An = [pDwag1 — a,y?)/2mekT?|w,? + 


where ky = 0 for a = Oand 1. In order to determine the j 
fluence of the reaction 4, — 24 on the transport coefticients , 
A, Dy», numerical calculations are performed 
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Comments on ‘‘Characteristics and Sound 
Speed in Nonisentropic Gas Flows With 
Nonequilibrium Thermodynamic States” 


Harold Mirels 
Lewis Flight Propulsion Laboratory, NACA, Cleveland, Ohio 
February 24, 1958 


ao 1 FOLLOWS the viewpoint of reference 2 and de- 
fines the local ‘‘speed of sound” in any inviscid flow process 
by the relation 


ay? = (Dp/Dt)/(Dp/Dt (1 
where D(__)/Dtis the substantial derivative 
)/Dt = )/dt] + )/dx;] 


Eq. (1) reduces to the usual definition a? = (dp/dp) isentrope for 
flows wherein each fluid particle maintains a constant entropy 
This equation is intended to be applicable to nonisentropic flows 
such as those with external heat addition? and nonequilibrium 
thermodynamic states.!_ Reference I uses a» to establish the char- 
acteristic directions—i.e., Mach lines—at each point in the flow 
and proposes that the numerical integration of the equations of 
motion can be based on a characteristics net employing these direc- 
tions. The purpose of the present note is to question the general 
validity—i.e., usefulness—of Eq. (1) by showing that Eq. (1 
does not give the correct characteristic directions for two-dimen- 
sional flows with external heat addition. 

Consider the problem of the steady two-dimensional flow of an 
inviscid nonconducting fluid with external heat addition de- 
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fined by Q Qix1, x2). A perfect gas is assumed so that the 
necific heats ¢ and c, are constant and the equation of state 


i isp = pRI The equations of motion can be written 
+ p(Ouj;/Ox;) = O (2a 
puj(Ou;/Ox;) + 6;;(OP/Ox;) = O (2b) 
— | 1p) (O0p/Oxj) — (1/p) (0p/Ox;)] + 
p(ou;/Ox;) = pQ (2c) 
where? = 1, 2,7 = 1, 2. and repeated indices are summed. Here 


Ohas the units of heat energy added per unit mass per unit time. 
Methods for determining the characteristic directions of Eqs 
2) are given in various texts. It can be shown that only the 
coeflicients of the partial derivatives enter iato the determination 
of the characteristic directions (31.02) 
Hence pQ is not involved in this de- 


See, for example, Eqs 
and (31.05) of reference 3 
termination and the local characteristic directions are the same 
as those for isentropic flow—i.e., Q = 0. The local characteristic 
directions are thus defined by the ‘conventional speed of sound” 


for the class of flows described by Eqs. (2) 


Now, let us determine ay for this class of flows. Eliminating 


du; dx; from Eqs. (2a) and (2c) yields 


uj(Op/Ox;)] = }1 + 
[(Q/epT)/(u, 


2 [uj(Op OX 
p)(Op/Ox;)}} (4) 


Clearly a) does not identically equal a except for the special case 
Q = 0. Hence the characteristic directions based on ay are not 
the correct ones for the class of flows described by Eqs. (2) and 


the viewpoint of reference 1 must be re-examined. 
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Author’s Reply 


E. L. Resler, Jr. 


Associate Professor of Aeronautical Engineering, 
Cornell University, Ithaca, N.Y 


April 6, 1958 


D* MiIRELS’ NOTE takes exception to defining sound speed 
generally by ao? = (Dp/Dt)/(Dp/Dt) as in reference 1 
To make 
interpretation easier, first determine Q in terms of more meaning- 
Multiply Mirels’ Eq. (2b) by u; and add 
and thus obtain 


Allow me to interpret Mirels’ case more carefully. 


ful flow parameters. 
it to Eq. (2e¢ 
u;(0/Ox;) + (1/2)u,2 + (p/p)] = Q (1) 

If C7; + (1/2)u;2 + p/p = c,T,, where T, is the usual stag- 
nation temperature, we can replace Q in Mirels’ equations by 
u;(0/Ox;)c,T 
tions are not determined as easily as Mirels believes when he re- 


Written in this form the characteristic direc- 


tains only Q and claims Q does not involve any derivatives 


= 4 Op Oxj Op/ Ox 


Using the more general definition ay? u; ; 
one indeed obtains Mirels’ Eq. (4), and this becomes, with Q = 


ao? = a2} 1 + [(uj;/T) (OT,/dx;)/(uj/p) (Op/2x;)]{ (2 


Written in this form, it is clear that the value of aj) depends on 
whether or not 7, changes while the fluid particles pass through 
the wave front 


Mirels’ result implies that the sound wave is so 
thin that 7, does not change during the time it takes the fluid 
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particles to transverse the structure of the wave. With this inter- 
pretation, my definition? of sound speed agrees with what Mirels 
believes to be the correct expression. However, my definition 
also includes the more general case where heat addition does 
change the stagnation temperature of the fluid particles traversing 


a sound wave 


REFERENCE 


1 Resler, E. L., Jr., Characteristics and Sound Speed in Nonisentropic Gas 
Flows with Nonequilibrium Thermodynamic States, Journal of the Aeronautical 


Sciences, Vol. 24, No. 11, pp. 785-790, November, 1957 


The Injection of Air Into the Dissociated 
Hypersonic Laminar Boundary Layer 


Sinclaire M. Scala 
Research Engineer, Missile and Ordnance Systems Department, 
General Electric Company, Philadelphia, Pa. 


March 7, 1958 


SYMBOLS 


©4 = mass fraction of atoms 
fy = pwldue/ dx) | similarity stream function) 
h, = enthalpy of the ith species, including chemical 
mu Putw Mass-transfer rate 
u = components of velocity, parallel and normal to the surface 
z absolute velocity of the ith species 
ut 4 = net mass rate of production of atomic species per unit volume by 
chemical reaction 
w(Ou Ov) shear stress 
= density 
“ = viscosity 
Subscripts 
e = edge of boundary layer 
u = wall 
” = denotes differentiation with respect to » 


T FIRST APPROXIMATION, dissociated air may be treated as a 
binary mixture of air atoms and air molecules.! In order 
to include the effects of mass transfer into the boundary layer, 


Pe = 5.7 ATM Tw = 3500°R 


-------- EQUILIBRIUM 
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Fic. 1. Boundary-layer profiles for the frozen and equilibrium 
boundary layers (altitude 100,000 ft., Mach Number 20). 
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it becomes necessary to introduce a third chemical species and 
hence a diffusion equation. We have avoided this 
complexity by considering the injection of air molecules into the 
boundary layer, and hence the theoretical treatment is accom- 
plished within the framework of a binary mixture gas. 

The stagnation-point boundary-layer equations constitute a 
coupled seventh-order system of nonlinear differential equations,? 
and have been solved for variable Prandtl and Lewis Numbers 
by means of an IBM 704 digital computer for two limiting 
chemical constraints—frozen flow (wa = 
catalytic wall, “‘local-equilibrium” 
Fig. 1 are shown the velocity, temperature, and atom-concen- 
tration for different air-injection rates. As 
expected, the major effect of injection is to thicken the boundary 
layer which reduces the velocity, temperature, and concentration 
gradients at the surface, consequently reducing the heat transfer 
Q and surface shear stress. 


second 


0) over a completely 


and flow (wa # 0). In 


profiles several 


With regard to the local-equilibrium constraint, it is noted 
that since the composition is constrained to be in equilibrium, 
the chemical-source term wa is not zero and was obtained by 
solving the diffusion equation for the atomic species. In this 
case, the effect of mass transfer is to alter the chemical-source 
a < ©) 
0) at a greater distance 
Thus, the 
bulge in the temperature profile, due to the local chemical heat 
release in the equilibrium boundary layer, 


distribution so as to cause the atoms to recombine (w 
and the molecules to dissociate (wa > 
from the wall than when mass transfer is absent. 


is not particularly 
diminished by the presence of mass transfer but merely occurs 


0- 
ALTITUDE: 100,000 FT 
MACH NO: 20 
2 -20+ = 3500°R 
fw=-2 
= 
-40} 
-6 + + + + + 
a5 10 15 20 25 


Variation of source distribution with mass transfer for 
the equilibrium boundary layer. 


Fic. 2. 


8.0 EQUILIBRIUM FROZEN (Wa = 0) 


| | ALTITUDE: 100,000 FT 
MACH NO. 20 
7.0 T 
N 
6.0 
5.0 
>. N 800 R 
oO 
3500°R'N 
Tr 
zie 3.0 
2.0 
0 a 4 6 8 1.0 1.2 1.4 
s 
x10" FT? SEC 


[due 
ax 
Hypersonic laminar boundary-layer heat transfer vs. 
mass transfer. 


Fic. 3. 


Sr ACE SCIENCES—JULY, 1958 
140 7 
_| 
ALTITUDE 100,000 FT 
MACH NO. 20 
ion 
o 
2 Saf 1,,-3500°R 
x 60 — 
SNR 
— FROZEN ( W,=0) 
.2C | 
10) a 4 6 10 12 14 
Mw Fs 
* O sec SEC. 
Fic. 4. Hypersonic laminar boundary-layer shear function y 


mass transfer. 


at a greater distance from the surface. 
the temperature gradient at the 
the heat transfer. 
In the absence of mass transfer or slip at a surface, the energ 
flux Q with respect to the mass average velocity is, of cours 
> 


The net effect 
surface 


is to reduc 


and, consequent); 


identical with the heat transfer to the surface since the maer i 
scopic velocity is zero at the wall. If mass transfer | 
then the energy which passes through the wall at 
stagnation point is given by (neglecting radiation): 


Qu = [K(OT/dy) — pi vi hil (1) | 


stream 
occurs, 


and includes the contributions of conduction, diffusion, and 
fluid injection. The latter expression was normalized by divid- 
ing by the square root of the stagnation-point velocity gradient } 
and is shown plotted in Fig. 3. 

In addition, the effect of mass transfer on the skin friction is | 
shown in Fig. 4, where 7, is the surface shear stress. ; 


As in the case of no mass transfer,? when the wall is full 
catalytic, the effects seem to be independent of the rate of ; 
homogeneous atom recombination since the percentagewist 
reductions in heat transfer apply equally well to either of the 
limiting chemical constraints which presumably bracket th 
actual nonequilibrium state. 
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Flow Against a Vertical Plate With Large 

Suction 
J. R. Foote 
Division of Engineering Sciences, Purdue University, 

West Lafayette, Ind. } 
March 7, 1958 


N A RECENT NOTE! the author presented a solution of this prob- 

lem for the case of small values of a suction parameter. That 
solution appears to be exact, given enough terms of the series, 

but this is not proved. However, early truncation of the series 
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forced the use of only small values of the suction parameter ¢ 


Symbols defined in the previous note.') Furthermore, ¢ 
entered the series in such a way as not to permit an asymptotic 
f the series in powers of 1/ go. 


exp insion © 
while studying a perturbation technique with other 


However, 
nonlinear problems, a ‘solution’ was found which seems to hold 
for arbitrarily large go. Indeed, graphs of this function and its 
first derivative have just the behavior one would expect as ¢o 
grows 

The approximation to a solution valid near & = 0 is found by 
perturbation around a starting function which satisfies boundary 


conditions at € = 0. Into the mathematical problem 
Lig) + ee” — 1) 
gi) = go, ¢'(0) = (), g’(o) = 1) 
js inserted the function 
g(t) = + 


where //(£) is assumed to be known. The resulting linearization 
gives 

+ Hh" — 2H'h' + = —L(H) (2) 
Clearly the first approximation ts 


H(t) = go = constant 


The problem becomes 
+ gh" = -1 


h(o) = h'(0) = 0, ~) not specified 


Itisa matter of elementary integration to find 


o(t) = oo + — — (c/ go?) (1 — e~™) (3) 


where c should be selected to satisfy ¢’( 2) = 1. However, 
when ¢'(£) is examined, it is seen that such a ¢c cannot be selected 
to satisfy the condition exactly. One choice is better than all 
others—-¢ = g since this for < gy will give curves which 
satisfy any requirement for engineering data near — = 0. More 


precisely we should have &< VY 2g). For this choice of c, 


term begins to pull the curve away from its 


Ultimately the 
asymptote and, heuristically, one might merely delete the term 
at that point. This desirable deletion may be based more firmly 
by examining the formal solution valid near — = ©, as seen 
below. 

If a further effort to better Eq. 3 is made by choosing this 
function as a new //(£), a certain difficulty arises. The resulting 
*§ making analytical treat- 


Eq. 2 has coefficients containing e7 
ment very difficult. Furthermore, if Eq. 3 is simplified to a power 


series and only terms of least order are kept, 


H(t) = og + Be, 


h(é) = 


But it turns out that the solution for h(&) gives only the next 
power of € neglected in expanding Eq. 3. Thus, the idea of suc- 
cessive perturbations can succeed here only if the function e~ * 


can be retained in the coefficients of Eq. 2. 


The behavior of solutions to Eq. 1 near & = © is examined 
by using s = 1 /¢in Eq. 1 
” 
Wig + — + 
— 2g)¢’ + sig’? —1=0 
¢ =, lim [—s%¢’(s)] = 0, lim [—s?¢’(s)] = 1 
«x 


Linearization as before gives 


+ — + 26s — 2H + — 
3(2H’ +sH"')h = —M(A) (5) 
The condition at ¢ = 0 is satisfied by H = 1/s; also 


H( ~) = (0), 


h( ©) = (6) 


lim [—s*h’(z)] = 0, 


ot) 
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Eq. (5) gives 
sth’’’ + (62° — s)h"’ + (622 — 4h’ = 0 


Although this equation has an irregular singular point at > = 0, 
formal methods may be applied to find the general solution 


h(z) = + c[(1/3s3) + + [(1/323) 4 
(1/s) }@(1/sY 2) +6 [(1/62?) + (1/3)] — 
J Pi/sy/2 
where P(1/sY 2) = (2/y dt 
0 
The first condition of Eq. 6 forces c. = c, = 0, while the third 
condition gives ¢ = go. Thus, ¢(f) = & + go, is valid near — = 


This indicates again how to select ¢c in Eq. 3 and how to 
correct Eq. 4 after & + ¢ has been reached, which is assured for 
sufficiently large go, say go = 10. 

With suitable experimental apparatus, the value of each 
parameter entering the problem of viscous incompressible flow 
against porous cloth can be determined, and thus the theory 
the parameter 
It could 
be measured at a point Xx» on the x-axis, where the velocity varies 


tested. Experiment seems necessary to find a, 


relating to potential flow outside the boundary layer 
linearly; then a = — M/X». Equivalently, the pressure could be 
measured there and a Bernoulli equation used. Compared to 
kinematic viscosity v, the parameter a is likely to be large so the 
parameter a = Va ‘y will be large. Thus, x = £/a@ indicates 
that the already rather rapid change to potential flow in terms of 
& (Eq. 4) is even more rapid in terms of x. This will make detec- 
tion of boundary-layer effects near a porous wall difficult for large 
values of g. Furthermore, experiment to determine gp» seems 
necessary. One finds 
go = PV wxo/v 

Here P is 
used for relative porosity, essentially the ratio of speed down 
It is estimated that, 


and again viscosity magnifies the parameter size 


stream to the speed upstream of the wall 
depending on the cloth and air conditions, g» can range from 


zero to 10%. 
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Calculation of Axisymmetric Isentropic 
Spike Surfaces 


Ernest C. Kennedy 

Senior Research Engineer, Convair, A Division of 
General Dynamics Corp., Daingerfield, Tex 

March 7, 1958 


I’ A CONE is placed in a supersonic air stream at zero angle of 
attack, a conical shock wave will be formed at the apex under 
In the region back of the shock, the flow is 
isentropic and supersonic. The Mach lines in this region do not 
If the cone surface is altered, the flow field 
If the change in the surface 


certain conditions. 


converge at a point. 
about it will change according], 
direction is gradual, the flow may be isentropic, but the associated 
family of Mach lines (left characteristics) will not, in general, 
intersect at a common point. In designing a diffuser for a ram- 
jet engine it is often useful to deal with an isentropic flow field 
in which the Mach lines converge at a point. 

The object of this paper is to describe a method of constructing 
an axisymmetric spike surface which will foree the Mach lines 
in the characteristic region to coalesce at a specified point in such 
a way that the flow will be isentropic throughout. The problem 
naturally divides itself into several phases 

Since our surface is axisymmetric, it will suffice to consider a 
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Hence when we speak of the 
we mean the meridian of the spike, as shown in 


plane cross section of the spike. 
‘spike surface, 
Fig. 1. To start the solution, we take the origin of coordinates 
at the apex of the cone and choose a point P on the conical shock 
at a fixed distance (say, 10) from the cone axis. 
thru which the left running characteristics are to pass. 

From reference 2 we get the value of the Mach Number (.1/,,) 
and the flow inclination (6,) just back of the conical shock. The 
Mach angle (@,) is obtained from 1/,. Since the network is to be 
pushed toward the spike, it is obvious that the left characteristic 
thru P must start off with an inclination of Z, = 0, + a, + 
180°. 
so it is merely a case of finding the points of intersection of the 
first left characteristic with this series of rays. This is done by 
drawing the characteristic line in short segments with the end 
points on adjacent rays, taking care to make the inclination of 
each segment equal to the average of the values at the two end 


” 


This is the point 


The values of 1/7 and @ are known along each conical ray? 


This process is carried on down until the characteristic 
This line will then form the 


points. 
line meets the cone surface at Q. 
boundary between the conical flow field and the characteristic 
region. 

After determining the first left characteristic we return to the 
point P and impose a small incremental turning of flow in accord- 
ance with the Prandtl-Meyer expansion given below 


y= V4 + 1)/(y — 1) arc tan 


V ily — D/y + 1D) CAP — 1) are tan VM? 1 


where y is taken as 1.4. Since the parameters are known at P 
(for any given turning) and at a series of nearly equally spaced 
points on the boundary curve, we proceed to find the parameters 
at the intersection of the right characteristic (emanating from the 
boundary points with an inclination of Zr = 9r — apr) and the 
left characteristic thru ? by means of the standard network 


(1) 


formulas: 
M;* = — 0, + cot ag + cot az + [sin @sin a(d/r)|_, + 
[sin sin a(d/r) /} (cot a/M*), + (cot a M*)p} (2) 
6; = (cot M3* — cot a, + 
6, — |sin sin a(d/r) (3) 


where d = length of side of mesh, 7 = ordinate of midpoint of side 
of mesh, 1/* = 
the subscript 3 indicates the point to be calculated. 
Eq. (1) holds for two-dimensional, isentropic, irrotational flow of 
a perfect gas. 


Mach Number referred to sonic conditions, and 
Observe that 


TABLE 1 
Vp Zi = 6) + ap + 180° 
7.6880 2.6367 42.2341 22 .2880 209 . 9260 
8.8721 2.5816 41.0000 22.7899 211.6620 
5378 40.0000 23.2061 213.0782 


9.8721 2. 
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The second left runner is carried down until it intersects | 
surface. We then go back to P, impose another incremen 
turning, and start down the third left characteristic and proce 
In this way a series of points are located on the q 
sired isentropic spike surface. 


as before. 
This process can be carried 

until the desired degree of compression (measured by Ma 
Number) is obtained or until the surface attains the requir 
slope. 

The spike surface is made up of short segments drawn from oy 
left characteristic to the following one. For example, let 6, an 
7, be the last two points determined on the surface. The ney 
left characteristic is brought down until it intersects the lip 
thru 6, and 7, at the point 4. The value of 64 is found by inte 
polation between two points on the characteristic. A line | 
then drawn thru 7, with an inclination of (6;, + 94) /2, meeting 
the left characteristic at B. 
Three or four such iterations vield a ney 


Then 6% is averaged with 6, an 
the process repeated. 
surface point &,. 

As an illustrative example of turning, let us consider the cay 


where M/, = 3.0217, @. = 15°. At P we have 

A, = 2.6367 (reference 2) ; 

6, = 7.6380 (reference 2) 

= 42.2341 (Eq. (1) 

@, = are sin (1/.1/,) = 22.2880 (definition) 

Zi = 0, + a, + 180° = 209.9260 (definition) 

Now associate ¥, with 6,, letting ¥,, = 42.2341° be the refer 
ence 6. Let the turning be in increments of 1°, except in the first 
case where we turn 1.2341°. The results of the first few turnings | 


are shown in Table 1. 

The first row in Table 1 gives the initial characteristic or 
boundary curve. 

It is evident that a spike surface constructed as above will meet 
the required specifications. That is, if the spike surface is placed 
in a supersonic air stream at design Mach Number and a char- 
acteristic network is built on it, then the left characteristics will 
converge at a point and the flow field will be isentropic through- 
out the region back of the conical shock. 
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Flow Induced by Jets 


Sir Geoffrey Taylor 


Cavendish Laboratory, University of Cambridge, 
Cambridge, England 
March 7, 1958 


r A RECENT ARTICLE on flow visualization,! A. M. Lippisch 
shows the streamlines of the flow induced by a two-dimen- 
sional jet, and he remarks that when the jet emerges at right 
angles to a plane they appear to be parabolas. He compares 
them with streamlines of an irrotational flow 
known flow potential 7 = s!'? as the vortex sheet which repre- 
sent the jet boundary, we found that the theoretical computation 
of the outer flow is in close agreement with our experimental 
results. The use of the words ‘‘vortex sheet’’ here is perhaps 
a little misleading. The sudden change in flow velocity at the 
cut in the theoretical flow field which idealizes the action of the 
jet on the surrounding fluid seems to be more nearly related to a 
sheet of sinks than to a sheet of vortices. 
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TABLE | 


Stream Functions for the Flow Into Jets 


Jet Sink 

Tet Velocity Strength Two-Dimensional Stream Function or 
Type m Stokes Stream Function 
F2H = r'* [eos (1/2) @ — sin (1/2) 
F2l # ¥ = eos (1/2)0 
F3H Z const = cosé 
F3I Z const vs = r(l + cos 8) 
T2H const const v2 = recosé@ when 0 < @ 

» = —reosA when 0 — 
rt const const method fails 

P, 2(—cos 9) } 


During World War II, I had some connection with the proj- 
ect F.1.D.O. in which it was hoped to clear airfields of fog 
either by forced or thermal jets directed upwards. In that con- 
nection I made calculations similar to those of Dr. Lippisch and 
found, as he has done, that when a two-dimensional forced jet 
emerges from a horizontal plane the streamlines of the inflowing 
air consist of two sets of confocal parabolas with horizontal axes 
directed outwards from the jet, but when the jet is unbounded 
the inflow streamlines lie on one set of confocal parabolas with 
axis in the line of the jet but lying in the opposite direction. 

Besides the two cases considered by Dr. Lippisch, I also investi- 
gated the inflow into conical forced jets and thermal jets. These 
results were never published except in the form of internal re- 
ports? to the British Ministry of Supply. The appearance of 
Dr. Lippisch’s paper therefore seems a suitable opportunity to 
mention them. In Table 1 they are summarized and classified 
by means of six symbols, F or T meaning that the jet is forced or 
thermal, H or U meaning that it emerges from a horizontal plane 
or into unbounded space, and 2 meaning that the jet emerges 
from a slit or heated line or 3 when it emerges from a small round 


1-4 
1:0 4 
T3H 
r cos O 
heat source 
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-1:0 + T 
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+ 
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Fic. 1. Streamlines for inflow to jet rising from a heat source. 
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orifice or is formed at a heated spot. Dr. Lippisch’s two cases 
in this notation are F2H and F2U 

This classification covers eight cases 
using the idea that, so far as the inflow is concerned, the boundary 
of a jet may be considered as a surface over which air is en- 
trained at a rate proportional to the velocity of the jet. The 
breadth and velocity of the jet at any height was determined 
The jet 
was assumed very narrow, though that is not a necessary restric- 
tion, and the flow calculated was that due to a distribution of sinks 
over the plane of symmetry in the two-dimensional case or the 
axis in the three-dimensional case. In Table 1, the variation 
in jet velocity with height Z above the source is given in column 
2, and the variation in m, the strength of the sink distribution, 


All but one were treated 


through the equations of continuity and momentum. 


in column 38. The two-dimensional stream function ys or the 
three-dimensional Stokes stream function y¥, are given in column 
4, where ¢ refers to distance from the origin and @ angular dis- 
tance from a vertical plane or line passing through it, @ = 0 being 
directed vertically upwards 

The streamlines fof F2H are two sets of confocal parabolas 
with horizontal axes. The streamlines for F2U are a single set 
of confocal parabolas with common vertical axis pointing away 
from the jet direction. The streamlines for F3H are horizontal 
straight lines. The streamliaes for F3U are confocal parabolas 
derived from the velocity potential @ = sin? (1/2)0 

The jets produced by a horizontal line source of heat in a sur- 
rounding atmosphere of uniform potential temperature have 
the peculiar property that the vertical velocity ia them remains 
constant while their width increases uniformly with height. If, 
therefore, the jet starts from a horizontal plane, the inflow might 
be expected to consist of uniform currents entering with equal 
velocity from both sides. This flow Table 1 for 
T2H, but since the velocity does not vanish at infinity it may 


is eatered in 


not be the ouly possible flow 

The stream functions for T3H and T3U involve the use of the 
Legendre function of order (2/3). No tables of this seem to 
have been published, but the Cambridge Mathematical Labora- 
tory supplied me with values of P2/3 (cos @) aud (d/d0)P2/3( cos @) 
computed on their Edsac computer. Using these, the streamlines 
y = 1 for T3H and T3U were computed. 
Fig. 1. It will be noticed that T3H is the only one of the seven 
cases summarized in Table 1 for which the streamlines are con- 


They are shown in 


cave upwards in any part of their length 
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Some Experiments With a Resonance Tube ina 
Supersonic Wind Tunnel; 


Merwin Sibulkin* and Thomas Vrebalovich 

Jet Propulsion Laboratory, California Institute of Technology, 
Pasadena, Calif. 

March 11, 1958 


_ NOTE GIVES a brief description of some exploratory ex- 
periments conducted with a resonance tube. The air col- 
umn in certain configurations of a cylindrical total-head tube 
(resonance tube), with its downstream end closed, oscillates at 
large amplitude when it is placed in a high-velocity jet or wind 
tunnel. Sprenger! found that the base of a resonance tube gets 

+ This paper presents the results of one phase of research carried out at 
the Jet Propulsion Laboratory, California Institute of Technology, under 
joint sponsorship of the Department of the Army, Ordnance Corps (under 
Contract No. DA-04-495-Ord 18), and the Department of the Air Force 

* Now with Convair Scientific Research Laboratory, San Diego, Calif 
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Resonance tube with “ring airfoil” trip. 


Fic. 1. 


much hotter than the free-stream stagnation temperature of the 
jet of air which is impinging on it. Dailey? and Hartmann? found 
configurations which oscillated at high pressure amplitudes, but 
they did not look for or observe any overheating. 

A photograph of the forward portion of the resonance tube 
which was used in these experiments is shown in Fig. 1. The 
model consisted of a thin inner cylinder, 3/4 in. in diameter and 
10 in. long, which was mounted in an insulated housing and 
instrumented with thermocouples and a static-pressure orifice 
1/2 in. from the mouth of the tube. The base could be set at 
any position from the mouth of the tube toa depth of lO in. The 
base was either a microphone or a water-cooled diaphragm. When 
this model, with the ring-airfoil device removed, was placed in 
the 12-in. supersonic wind tunnel of the Jet Propulsion Labo- 
ratory at a Mach Number of 2.8, the air column resonated at low 
pressure amplitudes at the quarter-wavelength organ-pipe fre- 
quency of the tube, and the tube remained cool. However, when 
devices, which we shall call ‘‘trips,’’ were added to the nose of 
the model, the air column resonated at high pressure ampli- 
tudes, and the thermocouples near the base of the tube registered 
average temperatures which were well above the free-stream stag- 
nation temperature of the wind tunnel. The intensity of the 
pressure oscillation and temperature rise depended on the type 
of trip and, for a given trip, on trip position. A typical trip is 
shown mounted on the front of the resonance tube in Fig. 1, and 
all of the data presented in this note were obtained with this 
“ring airfoil” trip. 

The results presented in Fig. 2 show how the pressure ampli- 
tude at the base and the temperature distribution along the tube 
are dependent on the base position. 7) is the stagnation tem- 
perature; Py is the supply pressure; and p’ is the r.m.s. pressure 
on the microphone at the base. The resonant frequency with the 
base at the 10-in. station was 294 cps, which is about 10 per cent 
lower than the natural frequency at 7). Furthermore, both the 
maximum temperature and the maximum r.m.s. pressure occurred 
with a tube depth of 10 in. A plot of 7 — 7) along the tube for 
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Fic. 3. Temperature distribution for base 10 in. from entrance 


of tube. 


the base at the 10-in. station is shown in Fig. 3 for the following | 


2.8; (2) = 28 
but the microphone replaced by a water-cooled base; and (3) \/ = 


0.66 (obtained by operating the wind tunnel in an ‘‘unstarted’ 


cases: (1) test-section Mach Number J/ = 


condition). Case (2) indicates that the air column is heated all 
along the tube since cooling the base mainly affected the thermo- 
couples nearest it. Case (3) indicates that the flow need not be 
supersonic in order for the resonance phenomenon to cause over- 
heating. 

Oscilloscope traces of the instantaneous pressure at the micro- 
phone indicated a sharp pressure rise which suggests the presence 
of shock waves within the tube. In attempting to understand 
the phenomena reported here, interesting questions arise involving 
(1) the dissipative mechanism inside the tube, which might be 
shock waves, and (2) the coupling between the internal and ex- 
ternal flows. 
made by both of the authors. 


Further investigations of these questions are being 
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A Proof of Duhamel’s Analogy for Thermal 
Stresses* 


A. J. A. Morgan 
President, AER, Inc., Pasadena, Calif. 
March 11, 1958 


— THE COURSE of some work! on the numerical detet- 
mination of the thermal stresses within thick-walled solids of 
revolution (which will be reported upon more fully at a later 
date) a method of proof for Duhamel’s analogy was found which 
may prove to be of general interest. It is the purpose of this 
note to present such a proof and to indicate wherein it differs 
from those found in the literature. 


* The results presented herein constitute a small portion of the work per 
formed under a course of research sponsored by the Army Ballistic Missile 
Agency (ABMA), Huntsville, Ala 
necessarily represent those of ABMA. 

The author would like to acknowledge the help provided by Dr 
Fung (California Institute of Technology) who allowed him to use the ua 
published lecture notes on his course in Thermoelasticity and the encourage 
ment provided by Mr. E. A. H. Hellebrand (ABMA) throughout the work 
under the contract cited in reference | 
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T(xg,t) #0 


F(xa,t #0 


I 
I 4 


I 
Yap ~ Yap 


# 


Fic. 1. A graphic representation of Duhamel’s analogy 


Duhamel’s analogy can be of great importance in the solution 
of thermal-stress problems since, by use of the relations compris- 
ing it, it is possible to establish a correspondence between a 
thermal-stress problem for a given body and a fictitious stress 
problem for the same body at a uniform temperature. 

Duhamel’s analogy can be proved by working solely in terms 
{the displacements in a purely formal manner. Such a method 
f proof is well known (reference 2, p. 203) and it will not be con- 
sidered in this note. In another method of proof, one works di- 
rectly with the stress-strain relations and the equilibrium equa- 
Timoshenko (reference 2, p. 206), using this latter method, 
However, several of the extremely 
We shall draw 


tions 
has provided such a proof. 
crucial steps are omitted from his presentation. 
ittention to, and rectify, these deficiencies. 

For convenience in the following, all arbitrary bodies, for which 
corresponding quantities are distinguished by the superscripts 
| and II, will have exactly the same initial (unstrained) shapes. 
We shall prove that, for Duhamel’s analogy to hold, they must 
ilso have exactly the same final (strained) shapes. 

Consider a Body I, for which a thermal-stress problem is to be 


solved. Pictorially, it may be represented as shown in Fig. 1, 
where* 
T (xq, 0) = temperature distribution within Body I, not identically 


equal to zero 


o'a3ixa, t) = stress distribution within Body I, not identically equal to 


zero 

sla8lxe./) = strain distribution within Body I, not identically equal to 
zero 

Pa! va, t surface tractions on the surface S! of Body I, not identi- 
cally equal to zero 

Fa! (xa, t) body force distribution with Body I 

Xa, = (x, ¥,2,/) = rectangular-Cartesian space-time point within 
or on S! of Body I 

a, 3 = tensor indices ranging over 1, 2, and 3 


Note that we do not assume that the Pg! and F,' are vanishing 
quantities. Furthermore, we assume that the material compris- 
ing Body I is of the Duhamel-Neumann type; that is, it obeys 
the following strain-stress relations: 


yas = [(1 + v)/E]loag — — aT (1) 
or, the stress-strain relations: 
= [E/(1 + [E (1 — 2v)| X 
where 
8 = oup = oi + a2 + 093 = first stress invariant 
@ = linear thermal coefficient of expansion 
E = Young's modulus for the material 
» = Poisson's ratio for the material 
and a, E, and v are assumed to be nontemperature dependent. 
To prove Duhamel’s analogy, we proceed as follows. From 


Eq. (2) we note that if all the strains vanish, there will remain a 
field of compressive stresses equal to — agg’, see Eq. (3) below. 
Thus, the stress field ogg! + oa will vanish with the strains if 
the tensile stress field ogg’ is given by: 


* Note that, in rectangular Cartesian coordinates, no distinction need be 


made between covariant and contravariant tensors 


467 


FORUM 


= 1 — 2 V0, 

That is, we consider a Body II for which the stress field is 
= Caf’ + [aET' (1 — 2») (4) 


The problem for Body II is to determine the quantities F,'', P,'', 
and y@g in such a way that 7" can be made to vanish identically 

The equations of equilibrium for Body II, on assuming that the 
inertia terms can be neglected —i.e., quasi-steady formulation 


are: 


Oxg) + = Oxg) > 


[aE/(1 — 2v)}] (O07! /Oxg) + = 0 (5) 
But, for Body I we have 
das! Oxg + =) (6) 
therefore, Eqs. (4) and (5) imply that: 
F,"" = F,' — [aE/(1 — 2v)] (07"/dxq) (7) 


From the boundary conditions for Bodies I and II we can write 


Tas YB = ET (1 — 2») 

Gas ¥3 + — 2p) = + [@ET'/(1 — 2p) (8) 
I I 
since Pg = gag vg, and the vg are the components of the unit 


normal to the initial surfaces, S' and S"', of Bodies I and II. 
Now, we can correct the deficiency in the proof of reference 2 
in the following 


I 
by determining the quantities 7° and ‘eS 


manner. From the Duhamel-Neumann relations, Eq. (1), we can 
write: 


I ~ I pil is 
= {(1 + v) E|oas' E) —al 


Yas. T bas (9) 
Impose now the restriction that, up to a rigid-body motion, the 
strained shapes of Bodies I and II must be exactly the same 
This means that the yag'! must be identically equal to the yg’. 
Thus, Eq. (9) implies that the temperature 7" in Body II must 
vanish identically. This completes the proof of Duhamel’s 
analogy. 

In summary, the relations comprising Duhamel’s analogy can 
now be depicted as shown in Fig. 1. Then, Duhamel’s analogy, 
Eqs. (4), (7), and (8), makes it possible to convert the thermal- 
stress problem for Body I into a pure surface-traction problem, 
with prescribed body forces, for Body II. Once the ae are 


obtained, the ogg’ can be determined from Eq. (4) 
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Stationary Conditions for Problems Involving 
Time Associated With Vertical Rocket 
Trajectories* 


Angelo Miele 
Associate Professor of Aeronautical Engineering, 
Purdue University, Lafayette, Ind 


March 18, 1958 


ONSIDERABLE ATTENTION has been devoted in recent years 
C to the problem of the optimum thrust programing along 
vertical paths. Trajectories of minimum propellant consump- 
tion have been investigated in references 1 to 4 for the case where 


* This research was supported by the United States Air Force through 


the AFOSR of ARDC under Contract No. AF 18 (603)-69 
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no time constraint is imposed. More general problems, involv- 
ing time, have been dealt with in references 4 and 5. 

With the present paper, the investigation* ® is extended. By 
techniques parallel to those developed in reference 6 for prob- 
lems not involving time, two first integrals are detected for an 
isothermal atmosphere in connection with variational problems 
involving time. An interesting consequence of this is a simple 
approximate relationship between altitude and Mach Number 
which holds (with a precision of 3-4 per cent) for all drag func- 
tions, regardless of the end conditions and of the analytical con- 
figuration of the functional form to be extremized. 

A rocket-powered vehicle traveling along a vertical path is 
considered. Its dynamical behavior is represented by the fol- 
lowing set of equations 


= 0 (1) 

Js=V+t+e4+ V) — V.8]/m} =0 (3) 
J, =8B-#=0 (4) 

Js =Bmar — B — 7? = 0 (5; 


where m is the mass, /: the altitude, V the flight velocity, g the 
acceleration of gravity (assumed constant), D (h, V) the drag, 
and V, the equivalent exit velocity of the rocket engine (as- 
sumed constant). The dot sign denotes derivative with respect 
to time. The quantity 8, which indicates engine mass flow, is to 
be such that 0 S 8 S 8B,,,q, at all points of the flight path. Notice 
that the above inequality is embodied in Eqs. (4) and (5) where 
— denotes a real variable and 7 another real variable. 

Among all sets of functions, m(t), h(t), V(t), B(t), &(f, and 
n(t) satisfying Eqs. (1) to (5) and certain prescribed end-condi- 
tions, the Mayer problem seeks that particular set such that 
AG = G; — G; is a minimum* where G = G(m,h, V,t). After 
introducing a set of variable Lagrange multipliers A,(¢) . . . A;(¢) 
the Euler-Lagrange equations’ are written as 


Ai + (As/m?) [D — V.8] = 0 (6) 

he — (As/m) (OD/dh) = 0 (7) 

hs + Ao — (A3/m) (OD/OV) = 0 (8) 
A — As( Ve/m) + Ay — As = O (9) 
ME = 0 (10) 

Asn = 0 (11) 


admitting, therefore, the following first integral 
+ + Ast [(V.8 — D)/m] — g} = (12) 
where C; is a constant. 
As the Euler Eqs. (10) and (11) point out, the solution are 


is discontinuous,t being generally composed of subares of three 
different kinds: 


n=0, = Bmaz (b) 
Ay = 0, A; = O— B= variable (c) 


For the subares of type (c) the Lagrange multipliers can be 
integrated in a closed form yielding 


A = C2V. exp [(V + gt)/Ve] (13) 
Ae = C2[(D/V.) + (OD/OV)] exp [(V + gt)/V-] (14) 
As = Com exp [(V + gt)/Ve] (15) 


so that the first integral (12) becomes 


mg + D{l — (V/V.)] — V(OD/oV) + 
(C,/C:) exp {—[(V + gt)/V.]} = 0 (16) 


* The symbol ¢ denotes initial point; the symbol / refers to final point. 

+ An analogous result can be obtained through the device of parametric 
representation of the engine properties, used by the writer in references 4 
and 5, 
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An isothermal atmosphere is now assumed—i.e., an atinosphe; 
where the speed of sound (a) is a constant. The drag is coy: 
quently rewritten as: 


D = (¥/2) puSCpo( A? exp ( —+gh/a?) 7 


where S is a reference surface, po the sea-level pressure, 5 + 
ratio of specific heat at constant pressure to specific heat »} 
constant volume, a the speed of sound, and Cpy the zero-lift dr 
coefficient. The latter is assumed to depend on the Mach Nyy 
ber WZ = V/a only. For the above type of drag function, 4 
variable-thrust subare admits the further first integral 
do(a?/yV.) — = C3 1s 
where C; is a constant—i.e., 
mg — (a?/yV.) [((D/V.) + (OD/OV)] + 
(Cs/C2) exp } —[(V + = 0 (19 
Elimination of the exponential term from Eqs. (16) and (| 
leads to the following result for the optimizing condition 
=2mg/ypS = — Vv] 
= Coo.M*[1 + + (d log Coo/d log M)} (2) 
W(M) = (Cpo/y) (M/4M,) [2 + (W/M.) + 
(d log Cpo/d log (2 
C; = Ci (C3 = Ci) 


bo 


the ratio of equivalent exit velocity of the rocket engine t { 
atmospheric speed of sound. Eq. has the form = (C,, 


where p is the local static pressure at altitude h and \/, = J, 


and represents a family of variable-thrust subarcs in the (¢, Mf | 
plane. Notice that there is one variable-thrust subarc for eact 
value of the constant C;. The particular value of C, associate 
with a given variational problem depends on the form of the 6-/ 
function and of the boundary conditions of the problem. Th 
case Cy = 0 corresponds to problems not involving time —ie 
to problems such that no time condition is prescribed at one of 
the end-points and 0G/oft = 0. 
Elimination of the mass from Eqs. (16) and (19) leads to 


(ygh/a?) — (gt/Ve) = (M/A) + log — + CG (4 


where C; is a constant. The above exact relationship is satisfie 
by all drag functions, and, in particular, by the so-called linear 
drag law! and quadratic drag law.! It also holds when the zer 
lift drag coefficient is proportional to an arbitrary power of the 
Mach Number.* ® 

Consider now the two terms on the left of Eq. (24) and calev- 
late the quotient (p) of their differentials 


p = d(ygh/a*)/d(gt/V.) = (25 


Assuming J, = 8, the above ratio vields p = 22.4 at / = 2 and 
p = 112 at VJ = 10. From this consideration it follows that 
for the all-supersonic domain, Eq. (24) can be approximated as 


ygh/a? = (M/M.) + log [@ — ¥] + G 26 


thus supplying a simple formula for the rapid evaluation of the 


h( M) relationship, valid for all arbitrary functions Cp(.M). Under | 
the same conditions of validity for Eq. (26), the time Mach | 
Number distribution can be approximated as 


ygt/a = (1/1) log M + (1/417) log — + 
M 
f [log (® — W)] + (27 
Mr 
where C, is a constant and .J/z a reference Mach Number. 


In closing, the author remarks that the approximations repre- 
sented by Eqs. (26) and (27) become unnecessary for the case 


where the relationship between zero-lift drag coefficient and ) 


Mach Number has the form Cpy = y 1/’~*, where y and x are 
appropriate constants. For such a case, in fact, the two fune- 
tions h(.M) and ¢t(.1/) can be integrated in a closed form. This 
circumstance has been demonstrated by Tsien and Evans fot 
x = 1 and x = 2 in connection with problems not involving 
time, and by the writer for x = arbitrary in connection with 
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hose more general problems where a time constraint is imposed 


in the flight trajectory. 
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Inviscid Hypersonic Flow Near the Stagnation 
Point of Oblate Ellipsoidal Noses 


Marcel Vinokur 
Associate Research Scientist, Missile Systems Division, 
Lockheed Aircraft Corp., Sunnyvale, Calif. 


March 18, 1958 


OR HYPERSONIC free-stream Mach Numbers, the flow in the 
F stagnation region of a blunt body is approximately incom- 
pressible. Lighthill'! obtained the inviscid solution near the stag- 
nation point of a spherical nose by assuming a spherical shock 
wave and expressing the incompressible, rotational-flow equations 
in spherical coordinates. This note presents the analogous solu- 
tion for oblate ellipsoidal noses which are of interest since they 
include spherical and flat-nose bodies in a one-parameter family. 

The elliptical coordinate system used is related to the eylindri- 


cal coordinate system by the equations 


(1 + (1 — 7’), = cé& (1) 
where the surfaces & = constant are confocai oblate ellipsoids of 


revolution, and the orthogonal surfaces » = constant are one- 
sheeted hyperboloids of revolution. The distance between the 
two focal points is 2c, and the positive z-axis is directed upstream. 

In the stagnation region, the incompressible continuity equa- 
tion for axisvmmetric flow is satisfied by defining a stream func- 


tion V such that 


ul —(OW/ On) V (1 + (£2 + n°) 
= (OW/dE)/YW (1 — + 9?) 


where u, v are velocity components in the &, 9 directions. 
be shown from the equations of motion? that W satisfies the equa- 


It can 


tion 
/OE*) — 9?) + /(1 + = + (3) 


where f(W) is determined from the boundary conditions. 
If the shock wave is assumed to be the ellipsoid = &s, then the 
boundary conditions are 


= —RU(1 — n?) (1 + &s7)/21 


(4 
(OW = —Uits(1 — 9%) 

and 
(W) = —(1 — + Es? + QW/RU + Es) 2 (5) 


Where l’; is the free-stream velocity and & is the ratio of the free- 


stream density to the density behind the shock wave. Since the 


incompressible approximation is valid only near the axis » = 1, 
the equation can be simplified by restricting the solution to 
Under this simplification Eq. (8) becomes 


+ 
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Fic. 1. Nondimensional stagnation-point velocity gradient vs 


density ratio across the shock wave for different values of 6/a 
/dE*) /(1 — 92) + (1 + &) = 
—(] — k)2U =? kl + 2 
The solution of Eq. (6) which satisfies Eq. (4) is 


w= + &) (1 — X 


}A[e/(1 + &) + tan-' — + Es*) — tan”! 

—(1 — k)*és?[BE? + 4 log (1 + &) — 8S (1 + &) - 
— 4 log (1 + + 8/(1 + + (7) 
from which it follows that the body nose is the ellipsoid — = &x. 


The boundary conditions (4) vield the relations 


A = (1 — 


k = [Bs + 2B, + — 48\(B, — 2B, +B) 


where 
(1 + &s?)B,; = és? [tan~' és — tan — + Ep x 
(15 + l0Es? + 3és4) + Es[BEs? + + + 
— 4 log (1 + — 8/(1 + En?) + 4 log (1 + Ez?) 
Bz, = 15[tan-! &s — tan~! + &s/(1 + és*) — 
te/(1 + (1 + = By — + /Es 


A quantity important for heat-transfer calculations is the 
nondimensional stagnation-point velocity gradient, given by the 


expression 


—(a/U) (dv/ds)yar = [(1 — R)*Es*(15Ee + 10ER* + 
/15R(1 + Es*)? — 24 1 + 


where s is distance from the stagnation point measured along the 
body surface and a is the length of the major semiaxis of the 
body-nose ellipsoid. It is plotted in Fig. 1 as a function of & for 
values of the aspect ratio of the body-nose ellipsoid 6 a ranging 
The solution 


from 0 (flat-nose cylinder) to 1 (spherical nose). 


obtained by Probstein® for b/a = 0 is shown in dashed lines for 


comparison. 


1 | 
) and (] 4 
| 
) 
Cs (4 | 
+ /30R(1 + — | 
Es(1 + Es? k)/2 > (S) 
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A Simplified Universal Rule for Subcritical 
Drag of a Supersonic Diffuser 


Charles Lee Dailey 
Chief Engineer, Wiancko Aeronautics, A Division of 
Wiancko Engineering Co., Pasadena, Calif. 


March 19, 1958 


SYMBOLS 


a = pressure 

= density 

if = temperature 

V = velocity 

D = drag on spillage stream tube surface 
\f = Mach Number 

6 = deflection through shock wave 

n = control surface inclination 

So = free-stream tube area 

Se = cowl entrance area 

Yo = radius of So 

Ve = radius of 

6 = stream tube radius having deflection, 6 
Cp = spillage drag coefficient based on S, 
( )o = free-stream conditions 

() = average value 


= SUPERSONIC DIFFUSERS Operating supercritically and with 
unblocked inlets (with all strong shocks swallowed), the flow 
field of the inlet is supersonic. Diffuser drag can be calculated 
exactly for simple shapes and with close approximation for al- 
most any case. In subcritical operation, however, the entering 
flow is mixed in character, and theory does not provide an exact 
solution for any case. 

Several attempts have been made to devise suitable methods 
for estimating subcritical drag. All of these methods use an 
application of the momentum theorem to a simplified approxi- 
mation of the flow field. It is observed’ that the drag rise in sub- 
critical operation depends linearly on spillage and is practically 
the same for two-dimensional plane-symmetric and axially sym- 
metric inlets without central bodies. Results are presented at 
Mach Number 2.0 only. 

From available drag data,‘* 
drag increment is roughly insensitive to inlet configuration and 
angle of attack and varies in an essentially linear manner with 
flow spillage for any inlet. It is suggested, therefore, that the 
simple theoretical method of reference 3 be used to estimate the 
subcritical drag rise for any inlet at any angle of attack. The 


7 it appears that the subcritical 


fir’ 
4, 
/ 


Fic. 1. Subcritical flow diagram for axially symmetric inlet. 


© PITOT INLET 

5 
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Fic. 2. Subcritical drag increment for any inlet. 


estimated drag rise thus depends only on free-stream Mach 
Number and flow spillage. 

For the sake of completeness in the present discussion, the 
method of reference 3 is developed in Fig. 1 for an axially sym- 
metric inlet with an arbitrary form of central body. 

The point of downstream sonic velocity will not always lie on 
the single shock ascending from the intersection point with the 
multiple-shock system as assumed in reference 3. A _ general 
point is therefore selected for convenience having a deflection j 
and located at the distance ys from the axis. 

Application of the momentum theorem to the annular region 
between the entering stream tube and the stream tube making 
the deflection 6 results in the following expression for drag due to 
gage pressure acting on the entering stream tube: 

—D — xr(ys? — (P — P.) = 

Vurl(ys? — cosé — VJ, 
where the bar over the symbols refers to suitably averaged values 
over the surface connecting the cowl lip and the shock wave. 

Solving for the drag and expressing it as a coefficient in terms 
of cowl frontal area (cookie-cutter area) and free-stream dynamic 
pressure one obtains 
Cp = 2[1 — (V/V.) cos 5] [(y5/y-)? — (¥o/¥e)?] — 

(2/yMo?) ((P/Po) — 1] [(v8/y-)? — 1] (1 

From the geometry of the figure it is evident that 

So/Se = (yo/Ve)? (2 


and from mass continuity it follows that 


pVrlys? — y-?) [cos (n — 5)/cos n] = po Vor(ys? — yw?) (8 


The latter expression can be rearranged to give 
[cos (7 — 5)/cos n] (p/p9) (V7 /V_) — (So/S.) (4 
= ( 


Ye [cos (n — 6)/cos n| (p Po) ( 4 Vo) as 


After substitution of Eqs. (2) and (4) in Eq. (1), one obtains 
2[1 — (S./S-)] 


Cp = = x 
[cos (n — 6)/cos (f/p.) (V/V.) — 1 


V cos (n — 6) p l P (5 
po Vo 


1 — — cos 
cos y¥M,? \ 
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READERS’ 


In Eq. (5) the velocity, pressure, and density variables with a 
har over them represent area-averaged values over the surface 
| .onnecting the cowl lip to the shock wave. Since the shock-wave 
pattern varies with mass flow, these variables must be mass-flow 
jependent, in general, and can be determined accurately only 
from a schlieren picture of the subcritical shock pattern. How- 
ver, for the special case of a pitot-type entry (axially symmetric 
sith no central body ), the variables are essentially independent of 
mass flow, and the drag on the entering stream tube depends 
inearly on free-stream capture area, S,. In fact, for this case the 
| following simple relation exists: 


S| = f(AL) (6) 


Cp —(S. 
This function is plotted in Fig. 2 and is compared to experi- 
mental data taken from references 4-7. 
It appears that this simple theory gives a good estimate re- 
gardless of inlet type, and, further, that accuracy of the estimate 
is possibly better than experimental accuracy. 
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On Certain Similar Solutions to Unsteady 
Laminar Boundary-Layer Equations in Low- 
Speed Flow* 


Kweng-tzu Yang 
Assistant Professor, Department of Mechanical Engineering, 
University of Notre Dame, Notre Dame, Ind. 


March 26, 1958 


UNSTEADY laminar boundary-layer problems have 
received wide attention both in compressible-flow cases! ~* 
and incompressible-flow cases.‘~7 The purpose of this note is 
to present an extension of the similar solution recently obtained 
by the writer? to the case where properties are variable and the 
Wall is porous and to the case of axially symmetric flow. For 
the two-dimensional case with variable properties, the low-speed 
laminar boundary-layer equations are: 


+ u(Ou/dx) + v(Ou/dy)] = 
(Ou.,/OX) + + (0/Oy) [u(Ou/dy)] (1) 


(Op/Ot) + [O(pu)/dx] + [O(pv)/dy] = O (2) 
+ u(OT/dx) + v(OT/dy)] = (0/dy) [R(OT/dy)] (3) 
with boundary conditions 


v, (x,t), T = T, (constant) 
(x,t), T,,(constant) 


y=0 u=0 


@ 


where w and v are velocity components in the x- and y-directions, 
respectively; ¢ the time variable; J the temperature; and p, yu, 


* This work represents part of research supported by the National Science 


Foundation. 
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FORUM 


Cp, k are density, viscosity, specific heat, and thermal conduc- 
indicates the local free- 
When the x, y, ¢ system 


tivity, respectively. The subscript = 
stream value and w the wall condition 
in Eqs. (1), (2), and (3) is transformed to the x, }, ¢ system ac- 


cording to 


and the continuity equation is eliminated by introducing a 


stream function y such that 


pu = 
py = —p,,[(Oy/Ox) + (O) /dx) + (O} /Ol 


together with Chapman and Rubesin’s viscosity-temperature 
law and assumptions that both c, and Prandtl number o are con- 
may be reduced 


stant, the boundary-layer Eqs. (1), (2), and (3 


to 


(O?y/OYOt) + (0?y/ox0 VY) — = 
[(po/p) (Ou, + + (4 
(00/dt) + (Q¥/dY) (00/dx) — (dy /dx) (00/91) = 
(Cv,./a) (070/O0¥?) (5) 
where v is the kinematic viscosity and 
(T — T,)/(T. — Ta), 
C= + 216)/(T. + 216 


In utilizing a second transformation of coordinate systems from 


(x, Y,t) to (x, n, where n = V(u./v,x) and assuming that 


f=y V =f(n), 
the following ordinary differential equations are obtained: 
—f +f" + 


+ (@/C) [f — (a@/2)nle’ = 0 (7) 


p) = al(n/2)f"’ — (p../p)] 6) 


under the specific free-stream condition 


= (x/L)/[1 — 


where u% is a constant reference velocity, L a characteristic length 
of the cylinder, and @ a constant flow parameter related to the 
unsteadiness in the undisturbed flow. The corresponding bound- 
ary conditions are 


=a, 


where a is a constant suction or blowing parameter. 
the condition f(0) = a implies that the suction or blowing velocity 
is a specific pure function of time. The density ratio in Eq. (6) 
is related to 6 by p/pj = 1+ |[(7./T~) — 10. It is noted that 
when property values are considered constant and a = 0, Eqs 
(9) reduce identically to those already treated by the 


)}=1, 00) =1, (9) 


f'@) = 0, f'(« 


Physically, 


(6), (7), 
writer.’ 

The analysis of an axially symmetric laminar boundary-layer 
problem in steady flow is generally of no particular difficulty 
since it may always be reduced to a corresponding two-dimen- 
sional problem by using the well-known Mangler transformation 
Unfortunately, this transformation fails to work in unsteady-flow 
Hence every axially symmetrical unsteady problem has 
Now if we consider an incompressible 


cases. 
to be treated separately 
flow over a blunt-nosed axially symmetrical body with constant 
property values, the momentum and continuity equations be- 


come, respectively, 


(Ou/Ot) + u(Ou/Ox) + v(Ou/OyY) = 


(Ou. + + (10) 
[O(ru)/Ox] + [d(rv)/dy] = 0 (11) 
where r = r(x) is the radius of the body of revolution. When 


the free-stream velocity distribution is again described by Eq. 
(8), which could be realized in the immediate vicinity of the for- 
ward stagnation point of the body, Eqs. (10) and (11) may also 
be reduced to an ordinary differential equation for a similar di- 


mensionless stream function F: 


F'" + 2FF" — +1 = 


al(n 2)F’ + 1] (12) 


= 
4 
Py 
Y= f, (p/p.)dy, 
| a 
= 
q 
(2 
) 3 
| 
(4 
| 
(5 
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where f = xV iu and the stream function W is introduced 
to eliminate Eq. (11). The pertinent boundary conditions are 
F(O) = a, F'(O) = Oand F’(~) = 1. It js interesting to find 
that when a = 0), Eq. (12) reduces to the well-known equation 
for steady flow. This analysis may readily be extended to the 
corresponding thermal boundary-layer equation, and also to the 
case of large temperature differences. 

As all these ordinary differential equations involve a number 
of parameters which could all vary over wide ranges, numerical 
solutions are now being programed on a differential analyzer, 
and results and their generalizations will be reported at a future 


date. 
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On ‘‘Buckling of a Simply Supported Plate 
Under Compression Reacted by Shear’’ 


Paul Seide 

Member of the Technical Staff, Space Technology Laboratory, 
A Division of The Ramo- Wooldridge Corporation, 
Los Angeles, Calif. 

March 28, 1958 


T A RECENT NOTE,! an approximate solution of the buckling of 
i simply supported rectangular plate under compression 
reacted by shear was presented. The writer would like to call 
to the author’s attention an earlier paper? in which a similar 
problem was treated somewhat more accurately and more gen- 
erally. In reference 2, the load system consists of unequal uni- 
form compressive forces applied to two edges of the plate and 
uniform reacting shear forces on the other two edges. In addi- 
tion, there are linearly varying shear forces along the compressed 
edges. The internal-stress distribution is then, exactly, for the 
problem similar to that of reference 1: 


VN, =0 
Ne y/a) 
Nyy = Ne ((b/2) — x] /a 


(1) 


\ 


The load system of reference 1 apparently differs from this in 
that there are no shear forces along the compressed edges, which 
leads most likely to infinite stresses at the corners of the plate. 
In the interior of the plate, away from the compressed edges, 
the stress distribution should tend to that given by Eq. 1, espe- 
cially if the plate is long. 

In his Rayleigh-Ritz analysis, the author neglects the shear 


JULY, 


forces in the plate interior and is thus inaccurate in his formyl 
tion of the analysis. A comparison of his results with thoge, 
reference 2 indicates that the conjecture of the last paragraph, 
reference 1, to the effect that some of the difference betweeg 
the results of the Rayleigh-Ritz analysis and a finite-differegg 
analysis might be attributable to this erroneous stress distributiog 
is correct. For the one case computed by both the Rayleigh 
Ritz method and the finite-difference method, the result of refer 
ence 2 falls between the author’s results. 

It further appears from the results of reference 2 that thg 
extrapolated portions of the curve given in reference | are erg 
neous. The portion tending toward a buckling coefficient off 
as the plate aspect ratio tends to zero approaches that value tog 
rapidly. This error is discernible from the curve itself since the 
calculated points indicate slower convergence to the valueg 
At the other end of the seale where the plate width becomes jg 
finite, since the buckling coefficient is expressed in terms of th 
width, it should become infinite rather than approach the valgg 
for pure bending. If, on the other hand, we interpret the caseg 
infinite aspect ratio as the plate length becoming very small, qe 
should still arrive at an infinite buckling coefficient since it be 
comes increasingly more difficult to buckle the plate. This com 
clusion is also borne out by the results of reference 2 for this 


special case. 
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Errata—‘‘Exact and Approximate Solutions for 
the Supersonic Delta Wing’”' 


L. R. Fowell 
Project Aerodynamics Engineer, Convair, A Division of General 
Dynamics Corporation, Fort Worth, Tex. 


March 26, 1958 


F D. HAINS of the Bell Aircraft Corporation has_ kindly 
e brought to my attention the omission of three symbol deft 
nitions in reference 1 (following Eq. 4.3.8) and reference 2 (fok 
lowing Eq. 4.4.26). These are 

3X2X3[.\ af*axi0X; 


A 


N = + 
Furthermore, a denominator term was omitted fiom the equa 
tion following Eq. (4.2.1) in reference 1 (Eq. 4.2.3 in reference 3) 
which should read 
= —[1 + (X2/tan 7)] sin 7 tan 8, 


The numerical results presented are not affected. 
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